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ANALYSIS OF CHARRING ABLATION WITH DESCRIPTION 


OF ASSOCIATED COMPUTING PROGRAM 
Fred W. Matting 
Ames Research Center 


SUMMARY 


A general method is presented for solving the problem of heat-shield 
response in the stagnation region of a charring type ablator. The analysis is 
actually for the stagnation point of an axisymmetric blunt body, but it is a 
valid approximate method for calculations in the stagnation region of any 
arbitrary blunt body. The analysis is applicable to wind-tunnel or flight con- 
ditions, and the heat loadings are either arbitrarily assigned or they are cal- 
culated concurrently with the heat-shield response. Surface heating (or 
cooling) mechanisms accounted for are those due to convection, radiation, 
homogeneous combustion, heterogeneous combustion, surface material removal by 
means other than combustion (includes erosion), and sublimation. Physical and 
thermodynamic properties of the ablating material are arbitrarily assigned so 
that calculations can be made for various materials. 

A typical application of the analysis is given as an illustration. The 
analysis is machine programmed for numerical solutions using a finite differ- 
ence scheme, and a family of computing programs is used. These programs are 
described and instructions are provided for using them. The programs can be 
obtained from COSMIC, University of Georgia, Athens, Georgia, 30601. 


INTRODUCTION 


Heat shields of the charring ablator type are in general use for frontal 
or stagnation region protection of current space vehicles. The charring abla- 
tor system has proven to be quite successful, at least up to entry speeds cor- 
responding to lunar returns. A considerable amount of performance data on 
various heat-shield materials has been and is being obtained in the laboratory, 
particularly in arc-jet wind tunnels. However, laboratory experiments cannot 
duplicate all the changing conditions encountered in entry flights. Thus 
analytical methods are needed for predicting heat-shield performance, and to 
evaluate the accuracy of the predictions, calculated results must be compared 
with laboratory data for any given heat-shield material. There are approxima- 
tions in the analysis presented in this report, and there are some uncertain- 
ties in the effective physical and thermodynamic properties of materials 
considered. (As an example, the thermal conductivity of some chars, which is 
temperature dependent, may exhibit hysteresis depending on the maximum temper- 
ature to which the char has been exposed.) When satisfactory checks between 
calculations and a range of experiments are obtained, we can assume that the 



material properties are being adequately represented, and we have confidence 
in making calculations for flights. 

A general method is presented here for determining the material response 
of a charring ablator in the stagnation region of a blunt body. The analysis 
has been formed for an axisymmetric stagnation point, but it is essentially 
valid in the stagnation region of any blunt body. The analysis is applicable 
to both wind tunnel and flight, and in either case the heat loading can be cal 
culated concurrently with the heat-shield response or arbitrarily assigned. 

The heat shield can have an assigned finite depth, or a calculation using a 
semi-infinite depth can be made. Surface energy transfer mechanisms accounted 
for include those due to convection, radiation, homogeneous combustion, heter- 
ogeneous combustion, surface material removal by means other than combustion 
(includes erosion) , and sublimation. Internally, the material receives (or 
loses) sensible heat, and it pyrolyzes (in a zone) from the original virgin 
plastic to (finally) pure char. 

An attempt has been made to represent the charring ablation problem 
mathematically in an accurate but simple manner. Some parts of the analysis 
are similar to the analysis of surface ablation in reference 1, but many of 
the basic mechanisms require a different analysis. The analysis has been 
machine programmed for numerical solutions. External conditions, material 
properties, and various calculation options can be read in arbitrarily, and a 
family of computing programs is used to handle the various types of cases. 

The results of a typical example are presented to illustrate the various cal- 
culations that can be made. The principal features of the computing programs 
are described, and instructions are provided for use of the programs. 


ANALYSIS AND METHOD OF SOLUTION 


Basic Approach and Approximations 

The analysis is concerned with the problem of charring type ablation that 
occurs in the stagnation region of a blunt body. The ablating material can be 
considered to be either very thick (semi-infinite) or of finite depth. The 
type of ablation studied includes pyrolysis or degradation inside the mate- 
rial. It also involves surface removal of material by heterogeneous combus- 
tion, by erosion or surface chemical reaction other than combustion, and by 
sublimation. Heat liberated by chemical reactions between the gases expelled 
and the external gas is also taken into account. 

The external heat loading (convective and radiative) on the ablator can 
be put into the calculations in two ways: either arbitrarily programmed as a 

function of time, or calculated concurrently with the heat-shield response. 

For wind-tunnel cases, the loadings, either programmed as input or calculated, 
are usually essentially constant. For flight cases, the loadings usually vary 
greatly with time. Their calculation involves the simultaneous solution of 
trajectory equations with the ablation equations. 
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The ablation response analysis is strictly set up for the stagnation 
point of an axi symmetric blunt body. The analysis is also used as a good 
approximation to calculate the heat-shield response in the stagnation region 
of an arbitrary blunt body. The external heat loading at the point for which 
the response calculations are made must be known (or calculated) . The analy- 
sis has also been used for response calculations at points far from the stag- 
nation region. This use is quite approximate, but it is not necessarily 
unreasonable when the heat loading at such points is relatively small. Since 
the analysis is essentially one-dimensional, the approximate use assumes that 
lateral heat transfer is negligible compared to heat transfer normal to the 
surface . 

As will be seen, the principal equation to be solved is a time-dependent 
energy equation that describes an energy balance in the interior of the mate- 
rial. A number of auxiliary equations are used to evaluate terms within it. 
The energy equation is a partial differential equation of parabolic type with 
independent variables, time, t, and depth in the material, y. Boundary condi- 
tions are determined by surface energy balance equations, and initial condi- 
tions must be supplied. The energy equation is solved numerically by finite 
differences . 


Conservation Equations 


The analysis is essentially a time-dependent energy balance along the 
stagnation center line of an axisymmetric blunt body. The one-dimensional 

coordinate system used is shown in 


- Char 


Pyrolysis 

zone 


4* Virgin plastic 



sketch (a) . The basic equations to be 
solved are the conservation equations 
for energy and mass, written in a sim- 
plified form that is valid in the stag- 
nation region. The densities of 
virgin material and pure char, respec- 
tively pp and p^, are considered 
constant. The conservation equations 
used are: 


Sketch (a).- Coordinates. 
Energy Equation 


where 


(ps c s + P g°pg) 
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+ Wgh sr 


v s “ v scomb + v sero + v ssub 

= v s (t) (i.e., independent of y) 
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( 2 ) 
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Wg rate of production of pyrolysis gases per unit space volume, g/cm 3 sec 

O 0) 

ihg rate of flow of pyrolysis gases past a point, taken as an absolute 
valued quantity, g/cm 2 sec 

h sr energy released by the pyrolysis reaction, cal/g of gas produced; 

h sr < 0 for endothermic reactions. This can be described as a func- 
tion of temperature, but it can often be conveniently approximated as 
a constant (see appendix A) . 

Equation (1) is considered to be a reasonable approximation. The thermal 
conductivity of gaseous material is neglected (or can be considered as lumped 
into the T( s term) . Also neglected are viscous dissipation and compression 
work. Secondary chemical reactions involving the pyrolysis gases as they flow 
through the solid material are not specifically treated. Secondary gaseous 
reactions, if known along with their rates, can be accounted for approximately 
by selection of the temperature enthalpy variation of the pyrolysis gases. 


Continuity Equation 


Pyrolysis Gases 

3pg 9mg 

9t 9y 


w„ 


Solid Material 


where 


3ps 9 f — ^ °Ps 
9t 9y ( p s v s) " ot w g 


w 


g = K (Ps ~ Pch) 


n 


The order of the pyrolysis reaction is n, and 




K = S r 


(3a) 


(3b) 

(4) 


( 5 ) 


is the Arrhenius formula for the reaction rate. The procedure will be to 
solve continuity equation (3b) for p s and substitute the solution into the 
energy equation (1) . 


Solution of the Continuity Equation 
We combine equations (3b) and (4) and have 

“fif = " K (Ps - Pch) n = of (p s " p ch3 (6) 

where Dp s /Dt is a substantial derivative (i.e., it follows a certain small 
mass of the solid material) . 
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Let 


J 


p s ~ p ch 
Pp “ Pch 


(7) 


where 1 - J is the degree of degradation. Then 


-(Pp - Pch) n ' llcjn CO) 

with the (usual) initial condition Ji(t^) =1. In the integration of 
equation (8), we define: 


z = J t K dt! (9) 

Equation (8) has the solutions: 

1 

J = £l + (n - 1) (p p - p c h) n_lz ] for n ^ 1 (10a) 

J = e“ z for n = 1 (10b) 

With the inversion of equation (7) , 

Ps = Pch + (Pp - Pch) J (y s *0 (11) 


The type of asymptote for J expected indicates that one should expect 
n > 1. In any case, in the computing program, J is not allowed to become 
negative . 

As noted above Dp s /Dt is a substantial derivative and it refers to a 
specific location in the material: 

y = y(.Yi> t) which is y = Yi + /f Vs dt l 

ti 

K = K(T) and T = T(y, t) 


K = K{T[y( yi , t), t]} 

‘i v 5 d tl , t)] 

This evaluation of K is used (in principle) in equation (9) , the integral 
for z. The solution for p s thus depends on the past history of a point in 
the material. When p s is determined, w g can be evaluated from equation (4) 
for substitution in the energy equation (1) . 

In the numerical computing program, z (eq. (9)) is obtained by a proce- 
dure that is equivalent to the equations developed above. Thus, 
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AX = X - X_! 

The subscript, -1, refers to the previous time line; and X and X_i represent 
total surface recession corresponding to the current and previous time lines, 
respectively. (Equation (59) is used to calculate X; see Surface Recession.) 
In the finite differencing scheme, the array, y(M), is used where M is the 
spatial grid point number in the spatial finite differencing. A new array is 
needed for following the time history of specific points in the material: 

y_l(M) = y(M) + AX 

where y^- 1 (M) for the previous time line and y(M) for the current time line 
represent the same point in the material. From the previous time line we have 
saved the arrays T_i(M) and z-i(M) (which are spaced according to the y(M) 
array). We interpolate against y_ i (M) to obtain the T_ 1 (M) and"z_i(M) 
arrays. We use equation (5) with the argument T_i(M) to obtain the i (- 1 (M) 
array. We then obtain 

z(M) = 0 . 5 [K (M) + KL i (M) ] At + zl i (M) (12) 

where At is the finite time increment between the previous and current time 
lines. This is the numerical solution to the integral in equation (9); the z 
values thus obtained are used in equation (10). (It will be seen later that 
T(M) and hence K(M) are known at this point in the calculation procedure.) 


Physical Properties of Solid Material 

Physical properties of the solid material will, in general, depend on 
temperature and on the extent of degradation. We need the specific heat and 
thermal conductivity of both the virgin plastic and pure char as functions of 
temperature. Data for partially degraded material would also be desirable, 
but in its absence we assume a linear relationship with density as below. 

c s = c ch + (Cp - c ch)^ P p _ ( 13a ) 

Thus : 

c s = c ch + ( c p - c chU (13b) 

Equation (13) is the general formula for the physical properties of the solid 
material. It is not perfectly consistent with a volumetric model sketched in 
the next section, but it is considered to be within the framework of 
approximations used. 

The formula assumed for thermal conductivity is somewhat exceptional; it 
is of the same form as equation (13) , but it is also modified by an empirical 
multiplying factor. 

% = [K ch + (Kp - K ch )J][l - 4(1 - Cj)J(l - J) (14) 
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The modifying multiplier in equation (14) changes the thermal conductivity in 
the pyrolysis zone (0 < J < 1) . The maximum modification (at J = 0.5) by Cj, 
a read-in quantity, is generally a reduction to account for the reduced con- 
ductivity in the pyrolysis zone of some materials. (Values of Cj of 0.6 to 
1.0 have been used in matching wind-tunnel data.) If no modification is 
wanted, Cj is assigned a value of unity. 

For some materials a reduction in thermal conductivity may be accompanied 
by a similar reduction in density in the pyrolysis zone, but since the overall 
effect on the calculations is small, it is neglected in the analysis. (The 
computer programs could be modified to include it.) Reduced thermal conduc- 
tivity across a narrow zone can have an effect because it influences the heat 
transfer across that zone. 

Physical properties are computed by a material properties subprogram for 
this purpose. Data may be either tabular or in equation form as inputs to the 
subprogram. 


Rate of Flow of Pyrolysis Gases 


The rate of flow of pyrolysis gases past a station, ihg, is a term in the 
energy equation (1) . This can be evaluated by summing up the rate of gas pro- 
duction between a station and the interior when one neglects the small change 
in gas stored in the interior. Equivalently, one can evaluate the rate of 
mass loss of solid material between a station and the interior to give the 
rate of gas production. This latter is the method that is used. The calcula- 
tion requires knowledge of the surface recession rate. This rate (the veloc- 
ity of solid material in the coordinate system. Vs) is obtained from 
equation (2) . The components in equation (2) , which are summed up, are 
obtained from surface calculations which are described below in the section. 
Surface Recession Velocity. 


The gas density, pg, in equation (1), is actually the weight of gas in a 
unit volume of space (with a portion of the unit volume taken up by solid 
material) . The gas is assumed to be at the local temperature of the solid 
material and at the external pressure, p s . Using the perfect gas law we have 


p g 


p s/g M g 

Rg T 


(15) 


where Rg is the universal gas constant and vg is the volume actually- 
occupied by the gas in a unit spatial volume. Using our simple model for the 
density of the solid material, we have 


Char, 

spatial volume 
fraction = l-J 



Virgin plastic, 
spatial volume 
fraction = J 


Sketch (b).- Volume fraction. 


p s = p ch (l - J) + Pp J (11) 

which can be visualized in a simplified 
way as shown in sketch (b) . We need 
the volume fraction actually occupied 
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where pg is the mass of gas per unit spatial volume. 

As noted above, the term riig in equation (1) is the mass rate of gas 
flow past a given station. When the small change in gas stored in a control 

volume is neglected, the gas flow rate 
past a station is approximately the 
gas production rate that occurs at all 
depths that are deeper in the material . 
In sketch (c) the dotted line repre- 
sents the control volume, and the 
solid and gaseous material is flowing 
from right to left. The term rhg is 
y BF a positive quantity. Continuity for 
the control volume with the semi- 
Sketch (c).- Control volume. infinite formulation is: 

AgCy, t) + V s (t)[p s (y MF , t) - p s (y, t)] = - g|- /^ MF P s d yi (I 9 ) 

Equation (19) , then, neglects the rate of change of the amount of gas stored 
in the control volume, which must be small. The value of y^p is constant as 
it is fixed in the spatial coordinate system. It is also large enough that no 
degradation reaches it, so it can represent infinite depth. 
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i — n— I'- stst- 


■s - r- 


1 I I 


I I I 



*g(y> t) = (p p - Pch^V^fJCy. U " O] - ^ / 7mf J(yi» t)dyij 

( 20 ) 

The finite depth formulation is similar but is slightly different in form. In 
place of equations (19) and (20) we have 

m g (y 3 t) - v (t)p s (y, t) = - J 7bf P s dyi (21a) 

y 

liigCy, t) = (p p - p ch )[v s (t)j(y J t) - T BF J(yi, t)d yi J (2ib) 

The limit, y^p, no1: in the coordinate system; it is the depth of the 

back surface of the material. It shrinks at the rate |v s |. Equations (20) 
and (21) are actually equivalent. As evaluated by the computing program, the 
integral in equation (20) or (21) is calculated for the current and the previ- 
ous time lines, and the partial derivative with time is approximated as a 
difference quotient with denominator. At. 


Surface Recession Velocity 

The surface recession velocity, v s , appears in equations (20) and (21) as 
well as in equation (1). For v ’ , we can rewrite equation (2) as 

= - ^7 O^scomb + ™sero + A ssub^ ^ 

where the m f s are > 0. With a charring ablator, there is essentially no sur- 
face recession prior to the beginning of pyrolysis. With the onset of pyrol- 
ysis (which proceeds back from the front surface) , a charred front face is 

quickly obtained. Thus, we write and use for v s : 

v s ^chcomb + ^chero + ™chsub) { 23 ) 


We now evaluate Achcomb’ ^chero* an< ^ ™chsub ■ (These are computed in the 
material properties subprogram.) 

Char com bustion rate- In calculating the char combustion rate, m c hcomb> we 
first evaluate the amount of oxygen that can diffuse to the char surface to 
perform combustion. In the model chosen, it is assumed that the pyrolysis 
gases burn, and have first claim on the available oxygen, leaving a lesser 
amount available for heterogeneous combustion. It is assumed that the homo- 
geneous burning occurs at a flame front near the outer edge of the boundary 
layer. The maximum air flow per unit area is 
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’“air " ( p air v air) max ~ (24a) 

max 

From boundary- layer theory we have 

’“air = ^ p air v air^e , = 

calc calc 

where <|> is a dimensionless stream function that normally takes the value, 
unity (see appendix B, eq. (B7)). For_ gases other than air, and for off- 
stagnation approximate calculations, 4> may be assigned values other than 
unity. The air available for oxygen diffusion is taken to be 



*air = f p air v air^ e = lesser of the values 

from equations (24) 


The available oxygen is 


’“ox ’“air^ox 


(25) 


(26) 


where C ox is the mass fraction of oxygen in the ambient gas . The oxygen 
theoretically used in the homogeneous combustion is: 

’“ox = ’“gw^ (27) 

theo 


where = ihg(0, t) is the rate of pyrolysis gas expulsion and E 2 is the 
stoichiometric ratio of oxygen to pyrolysis gas for the reaction assumed to 
take place. The excess oxygen is 


m, 


ox 

ex 


m. 


ox 


- m. 


‘ox 

theo 


(28) 


If m ox < 0, there is no excess oxygen and not all the pyrolysis gases are 
ex 

burned. The heat liberated by combustion of the pyrolysis gases is propor- 
tionately reduced, and according to our model, there is then no heterogeneous 
combustion. If 


then we use 


^ox 0 

(29a) 

ex 


m 0 x = 0 

(29b) 

ex 



in all subsequent calculations. For the concentration of oxygen by weight 
available for heterogeneous combustion, we can write 


J oxf 


j ox 



(30) 
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For heterogeneous combustion, oxygen must diffuse to the wall, giving a 
concentration at the wall, C QXW . The kinetic reaction rate depends on C QXW . 
Since there is no accumulation or depletion of oxygen in a quasi-steady state 
process, the diffusion rate and the kinetic rate of consumption of oxygen must 
be equal. 

For the oxygen diffusion rate, we use the Lewis analogy (with Le - 1) ; 
(see ref. 1). This states that the ratio of the rate of diffusion of a 
species to its concentration potential equals the ratio of the continuum con- 
vective heat-transfer rate (corrected by a blowing factor) to the enthalpy 
potential . 


^oxd 

^oxf “ ^oxw 


’K.c 

Ah 




( 


0.00836 



(31) 


where the continuum convective heat-transfer rate, q oc , the blowing factor, \p , 
and the enthalpy velocity V, are evaluated in the section. Surface Heat 
Transfer. With the evaluation for q oc (eq. (62)), and neglecting any vortic- 
ity correction on q oc , one can write 


_ ^oxd 
^oxf ~ C oxw 


0. 00836^4 



1.15 


*d 


(32) 


“oxd = K d ( C oxf ' C oxw) 


(33) 


where Kj depends on external conditions, nose radius, and the rate of expul- 
sion of pyrolysis gases. The theoretical maximum oxygen diffusion rate occurs 
when the wall concentration is zero. 


Aoxd - K d c oxf 
max 


(34) 


For the reaction rate of oxygen consumption in the surface combustion we 
use an Arrhenius type of rate constant. 

" T~ 


Kqx = K r e 


l w 


Then for a reaction of order j we have 


m. 


oxr 


" K r[ p air c oxwV 

L w J 


C E 

T w 


(35) 


(36) 


- v / Pt 2 C °X»Y 


T 

e w 


(37) 



where R a ^ r , the gas constant for air (or other external gas) , is approximated 
as independent of the oxygen concentration. Upon combining constants: 

K r 

Ki = - 2 - (38) 

d J . 

K air 


where is the constant read into the computing program. Then 


™oxr 




(39) 


The maximum theoretical reaction rate would occur with an oxygen mass fraction 
of unity, so we define 


K 


re 


= K 


(%)' 




(40) 


(Actually C oxw cannot exceed C ox f.) Then 

™oxr = ^re (^oxw^ 

With the quasi-steady state assumed, 

™ox = ™oxd = ™oxr 


(41) 


(42) 


With the substitution of equation (42) into equations (33) and (41), we have 
two equations in two unknowns, m QX and C QXW , that can be solved. In computa- 
tions, one generally assumes the combustion reaction to be of half order 
(j = 0.5). (This should be the case for the combustion of carbon with undis- 
sociated oxygen to carbon monoxide, C + (1/2)02 CO.) In two typical mate- 
rial properties subprograms, j = 0.5 is "built in"; in a third subprogram, j 
is an arbitrary read-in quantity (see appendix C) . Finally, the rate of oxy- 
gen consumption is converted into char combustion 


^ox^kdb 

m chcomb ” c x 


(43) 


where Cj^b is the reciprocal of the fraction of char that is carbon, and the 
noncombustible fraction of the char is assumed either to melt off or blow off. 
The quantity, C x , is the stoichiometric ratio of oxygen to carbon in the reac- 
tion. With the usual assumption of combustion of carbon to CO, C x is 4/3. 

For the theoretical maximum rate of char combustion by diffusion we have from 
equations (34) and (43) 

m oxd C kdb 

max .... 

m chcomb ~ r (44) 

d max x 
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From equations (40), (41), and (43), the maximum rate of char combustion, as 
determined by the reaction rate and the available oxygen, is 


K r e (C Q xf ) *^ C]<;db 

m chcomb “ 

r max x 


(45) 


The rate of char combustion, ^clicomb* cannot exceed either ra c hcomb or 

d max 

m c hcomb* q uantit y> ™chcomb* wou ld be the rate if the reaction rate were 

r max d max 

infinite, and ^ c hcomb wou ld be the rate if diffusion were instantaneous, 
r max 

These quantities are of interest, because the approach of ^chcomb to one 
value or the other denotes "diffusion control" or "reaction rate control" of 
the char combustion rate. 


Char erosion rate -For the char consumption rate by other than combustion 
or sublimation, we do not fully understand the mechanisms involved, and we use 
empirical calculations. We call this consumption rate erosion, 
although there may be chemical reactions (and energies) involved. The empiri- 
cal calculations should be related to wind-tunnel experiments in inert gases 
and also in air for a given material. The erosion rate will vary considerably 
for different materials because of the variation of density and structural 
integrity of the char produced. 

An equation that can be used for erosion is of the Arrhenius type: 

Ot o 

/Pt \ m - ^ 

m chero = p sw a lbf~) e w C 1 ~ K 2) (46a) 


Another type of equation that is based on experience at Ames Research Center 
has been used for several materials (e.g., for phenolic nylon): 


™chero Pswl- a l + a 2T w + a 3 C^\v 1 (1 ^ 2 ) (46b) 

Equation (46) (of whatever type) is in the material properties subprogram for 
a given material. The constants and m can be read in using open symbols 
available (see appendix C, Input Data), or the constants can be written numer- 
ically into the material properties subprogram. Threshold temperatures can 
also be used, above which the constants change values (so that wind-tunnel 
data can be fit very closely) . 

Experiments in inert gases and in air have indicated that for some mate- 
rials, the surface recession due to erosion and combustion may not be com- 
pletely additive (or the one may partly suppress the other) . The factor 
(1 - K 2 ) in equations (46) is used to account for this empirically, where 
0 < K 2 < 1 is a read-in quantity. As examples, for the Apollo type heat 
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shield material, K 2 = 0.6 has been used, while for high-density phenolic 
nylon, K 2 = 0 has been used. 

It is possible to write more complicated expressions for erosion than the 
examples in equation (46) by postulating several failure mechanisms. For exam- 
ple, a dependence on external pressure gradient would include a term that 

would be a function of p+- /R. 

2 

Sublimation rate- The calculation of the mass loss rate by sublimation, 
m c hsub> requires a bridging between the free-molecule rate and the diffusion 
controlled rate. This is developed in appendix C of reference 1. 

The free-molecule rate of sublimation in a vacuum is calculated by the 
Langmuir equation (refs . 1 and 2) : 


™chsFM ^cvPve^dy^kdb 




M, 


cv 


2TrR g T w 


(47) 


where C^y. is the pressure of a standard atmosphere in dynes/cm 2 , so that the 
vapor pressure, p ve , is in atmospheres. The constant, , converts the car- 

bon sublimation rate to that of the entire char. The assumption is that the 
carbon rate controls, and the other material in the char either melts off or 
is blown off. With evaluation of constants we have 


“chsFM 44 . 37A cv Pve^kd 


:b V T w 


The equilibrium vapor pressure, p ve , is evaluated as 


Pve = 


e 8 

e 



(48) 


(49) 


For the sublimation rate with diffusion control, we use the equilibrium 
vapor pressure, p ve , but consider that it can be modified to a value, p vm , by 
homogeneous chemical reactions. As in reference 1, this modification is 
approximated as 


* 2 . , 

Pvm \Pve/ 


(50) 


where E 7 = 1 represents no modification (which is the usual case) . We define 


P 




1 


(51) 


and we use the Lewis analogy equation as in reference 1 for the rate of mass 
diffusion (see also refs. 3 and 4). 
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(52) 


™chsd 


0.00836^q oc C kdb 

p © [v2 ‘ °- 00836 Vwi 


where converts from the carbon diffusion rate to the rate of entire 

char consumed, and where Mp is the molecular weight of the gaseous products 
(after heterogeneous combustion) . It is assumed that this is essentially the 
gas through which the carbon vapor must diffuse. 

For the mass loss rate by sublimation we have finally (see ref. 1) : 


1 

^chsub 


1 

^chsFM 


+ 


i 

"‘chsd 


(53) 


In the development above we have lumped all components of the char 
together, and we have assumed that the carbon sublimation rate controls the 
rate of consumption of char. The other components of the char may affect the 
constants used, such as equilibrium vapor pressure, p ve , and accommodation 
coefficient, A cv , but as a reasonable approximation, the constants for carbon 
are probably satisfactory for most chars. Typical values used for char con- 
stants are: Eg = 22.3, E^g = 91,600, A cv = 1.0, M cv = 33, and Mp = 15. Sub- 
limation, for many conditions, is a small factor in the ablation process. It 
can become important and even dominant, however, for very high speed flights 
such as from planetary returns. 


Quantities Required for Subsequent Calculations 


Total pressure- For the total pressure, p t ^ (atmospheres), for flight cal- 
culations we use a hypersonic approximation, twice the free-stream dynamic 
pressure with a correction factor (ref. 1) 


A X DV 2 
P t 2 = 101.3 


(54) 


where a typical value for Ai is 0.95, and the number, 101.3, accounts for the 
units in the equation. For the wind-tunnel case, p^^ is a read-in quantity. 
Also * 

p S2 = 1.013xl 0 6 pt2 (55) 

where p s ^ is the pressure in dynes/cm 2 . For the arbitrary heating rate com- 
puting programs, p t2 as a function of time is read into the program. For off- 
stagnation approximate calculations. A* should be adjusted so that p t in 
equation (54) becomes a local pressure, p w . This also applies to arbitrary 
heating rate calculations because, in these cases, equation (54) is inverted 
to obtain the free-stream density, D, from the programmed pressure, p t ^ or p w . 
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Enthalpy velocity -In place of using the enthalpy of the external gas for 
surface heat-transfer calculations, it is convenient to use an "enthalpy veloc- 
ity," in km/sec, which is defined as (ref. 1 ) : 

V 2 = 0.00836h st (56a) 

V = 0.0915/h^ (56b) 


where h st is in cal/g. 

Specific heat of exter nal gas-It is convenient to use an average specific 
heat for the external gas (as in ref. 1) . This is defined as 



T 

o 


c p ^T i 


We use the evaluation 


c p “ Eio + ^ 3 T w 


(57a) 


(57b) 


where Eiq and E 3 are normally constant read-in values. If a higher order 
representation is desired, Eiq and E 3 can be made functions of temperature in 
the material properties subprogram. Values of E 3 = 0.84xl0 _i+ and 
E 1 0 - 0.1267 have been used for air at flight conditions. 

Surface recessio n-The surface recession, x, is a quantity of interest, 
and it is also used in evaluating other quantities. 

X = |V S I (58) 


X = I 1 X dt 1 

t. 

1 


+ X, 


(59) 


Nose radius -The effective nose radius, R, is needed for the calculation 
of the convective continuum heat-transfer rate, q oc . The R is calculated as 
a quantity that has an empirical variation with surface recession (as in 
ref. 1 ). 


R = Ri [ 1 + E 17 X 2 + c x (e C 2 X - 1)] 


(60) 


The constants should be evaluated experimentally if possible; otherwise, a 
recession shape geometry must be assumed. For a negligible shape change, the 
constants have the value zero. In the arbitrary heating-rate computing pro- 
grams, the R, as a function of time, is programmed as an input. For approxi- 
mate calculations off the stagnation point, the value of R may be a 
fictitious effective value. 

Tumbling factor-Tumbling or oscillation in flight will reduce both the 
convective and radiative heat transfer into the point on the body that is 


16 


taken as the nominal stagnation point. The quantity, K tu , is used as a 
multiplying factor on the heat-transfer rates. It is evaluated empirically 
(as in ref. 1 ) as 


K 


tu 


K + (1 - K) 


1 



(61) 


where K is the fractional time that the point on the body (nominal stagna- 
tion point) is initially exposed to approximately stagnation conditions, and 
Xi and E 16 are values selected to give a realistic damping to the oscillation 
during an entry flight. For the entry of a small body that tumbles initially, 
but achieves stability, the values have been used: K = 0.25, Xi_j= R/10 cm, 

and Ei g = 1.0. For a nontumbling body, K is unity, giving a K tu of unity. 
For the arbitrary heating-rate computing programs, X tu is not used, but it is 
considered to be "built in" to the assigned heating rates. 


Surface Convective Heat-Transfer Rate 

The equations listed below are identical with equations (15) through (25) 
as developed in reference 1 . 

These equations are not used in the arbitrary heating-rate computing pro- 
grams; for those cases we use only equations (72) (inverted) and (73). 

In the calculated heating-rate computing programs, for laminar continuum 
flow, the surface convective heat transfer is evaluated as (refs. 1 and 5): 

q OC = A 4 J| V 1 - 15 ^ 2 - 0.00836c T](l + (62) 

1 K 1 \ yfDVR/ 

where A 4 is a constant, and Cg is a (generally small, 0 < Cg < 0.2) vortic- 
ity correction. (A 4 and Cg are not used with the arbitrary heating-rate com- 
puting programs.) For wind-tunnel tests, A 4 should be evaluated with a 
calorimeter; the value of A 4 for earth entries is approximately 1.1. With a 
blowing correction, we have 

V = ^oc C fi 3) 

where ip is evaluated in the section. Blowing Factor. 

For rarefied conditions, the evaluation of the surface convection rate in 

free-molecule flow is needed. We use the Newtonian type approximation (refs. 

1 and 6 ) : 

«FM * TOSS? t V2 - 0 -°0» 36 ¥«] < 64 > 
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The accommodation coefficient, A C q, is < 1 . If unknown, it is usually 
assigned the value unity (at or near the stagnation point) . In the transition 
regime between free-molecule and continuum flow, we use a relationship that 
gives a value of convective heat transfer that is bridged between the free- 
molecule and the continuum evaluations (ref. 1) . 





e 



(65) 


Equation (65) is used for all regimes; it automatically gives the correct 
free-molecule and continuum asymptotes. When the tumbling correction is 
applied, the convective heat transfer at the front surface is evaluated as 


% ^^ w K tu 


( 66 ) 


In equations (62) to (66) , all the convective heat-transfer rates that 
are actually needed are obtained. However, we can also calculate several 
additional related quantities. When there is no mass transfer (i|i = 1), 
equation (65) specializes to 


^oo 


V3C 



and equation (66) becomes 


% %o K tu 


(67) 


( 68 ) 


With mass transfer, we can define a modified \p that operates on q OQ 
(instead of on q QC as developed above) . The modified blowing factor is 
defined as (ref. 1): 


= 


^w 

%0 


(69a) 


^w 


= ^OO 


(69b) 


Multiplying equation (69b) by K tu and using equations (66) and (68) we have: 

q w = H 0 ( 70 0 

We substitute equations (63) , (65) , and (67) into (69a) to obtain an 
evaluation of \p. 
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e 


* 

? = 


1 



( 71 ) 


The development from equation (67) to (71) is an alternate form of convective 
heat-transfer analysis. It is entirely equivalent to the previous development 
in equations (65) and (66). The new quantities obtained, q 00 , and are 
considered to be of conceptual interest. Both developments are calculated and 
the quantities obtained are read out by the computing programs. 


Another quantity of interest that is calculated is the cold-wall convec- 
tive heating rate, q ocw (not in ref. 1) . This has been defined to include the 
oscillation correction, so it is the heating rate that would be received by a 
cold front surface with no mass transfer occurring. 


^ocw 


V 2 - 0.00836c p T w 


(72) 


For normal (nonrarefied) wind-tunnel conditions, the bridging relations 
between the free-molecule and continuum regimes are not needed; also there is 
no oscillation. We have then: 


K tU - 1 

(73a) 

¥ = Ilf 

(73b) 

% “ ^oc 

(73c) 

q w = = ^oc 

(73d) 


The cold-wall convective heating rate is given by equation (72) for this case 
also . 

When the arbitrary heating rate computing programs are used, the cold- 
wall convective heating rate, q QCW , as a function of time, is programmed as an 
input. The oscillation correction and rarefaction effects are already 
entered. We invert equation (72) to solve for q Q , and we use equation (73) 
to solve for the front-face convective heating rate, q w . 


Blowing Factor 

As noted in the preceding section, for the evaluation of the surface heat 
transfer by convection, it is necessary to know the blowing factor, \p . This 
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is evaluated empirically in terms of the blowing parameter, B, which is 
defined below. The total mass loss rate of material is 


™ T (t) = ihg(0, t) - v s (t)p s (0, t) (74a) 

P s (0, t) 

m T (t) = m g (0, t) + Oii chcomb + * chero + A chsub ) — — (74b) 

We also calculate, as a quantity of interest, the total mass loss of material 
(per unit area) as : 


Mp 



ihp dt ^ + Mp _ 


(74c) 


For the gaseous material that contributes to heat blockage one should add 
the oxygen that is consumed in heterogeneous combustion. The oxygen in the 
homogeneous combustion is not included because we are using a flame front 
model that considers the homogeneous combustion to occur at the outer edge of 
the boundary layer. For the blockage-contributing mass loss rate: 

1 p s (o, t ) 

A c = *g(0, t) + [A chcomb (l + C x ) + A 3 m chero + m chsub ] — 


(75) 

Division by the factor C^b means that only the carbon in the char contrib- 
utes to blockage. The material other than carbon is assumed to either melt 
off or be blown off as a solid. The constant, A 3 , allows the option of 
selecting whether or not the eroded carbon contributes to blockage. When 
A 3 = 1, the eroded carbon is included; when A 3 = 0, the eroded carbon is not 
included. The latter is probably the usual case, where eroded material is 
removed in the form of solid particles. A 3 = 0.5 has also been used which 
means that half the eroded material is assumed to be gaseous. A value of 
A 3 > 1 can be used to represent removal of material by some chemical reaction 
other than combustion and release of the products in gaseous form. The 
blowing parameter is defined as 

B = m c — (76a) 

^oc 


B 



1 

0 . OO836A4V 1 • 1 5 


In going from equation (76a) to (76b), the vorticity correction in 
equation (62) is neglected. 


(76b) 
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The blowing factor, ip, is approximated as a linear function of B 
(refs. 3, 4, 7) with a slope that depends on the molecular weight of the gas- 


eous products, Mp (after heterogeneous 
vented from going negative with large 



Sketch (d).- Blowing factor. 


combustion) . The value of \p is pre- 
B by arbitrarily connecting an 

"exponential tail" on the function at 
the arbitrary point, ip = 0.16. This 
has little effect on the heat trans- 
fer since \p is small over this por- 
tion of the range of large B 
values. An asymptote of 0.06 has 
been chosen for ip (as suggested in 
ref. 1) . The ip function is illus- 
trated in sketch (d) . The slope of 
the linear portion of the ^ curve, 
C^, depends on the read-in constant, 
A 5 , as well as on Mp . Values in the 
neighborhood of 1.7 have been used 


for A 5 (ref. 7); 
experiment. The 
power (ref. 7) in 
use : 


however, the slope is preferably evaluated by wind-tunnel 
Cw, dependence on Mp is taken as a negative one-third 
equation (77a) . T“ 


C,i, = 


For the ip function we use: 


S B 


other 

constants used depend on 

V We 

. A 5 



(77a) 

M p 7 3 



0.84 

s 



(77b) 

: 10 C^ 



(77c) 


for 

B s c * 2 

(78a) 

0.06 

for 

B > 

(78b) 


Radiation Energy Rate 

With the arbitrary heating-rate programs, the incoming radiation as a 
function of time is a programmed input. For the calculated heating-rate pro- 
grams, the gas cap radiation into the body is evaluated with an empirical 
approximation as in reference 1 (see also ref. 8 ) 

q R = E 4 RD E 5 V E 6 K tu (79) 

where the level of radiation and the surface reflectivity are included in E 4 . 
Typical values of the constants used for earth entries are: E 4 = 0.76*10~ 6 ; 

E 5 = 0.5; E 6 = 7.0. 
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For both types of computing programs, the net radiation into the body is 
the difference between q R and the reradiation out. 


^rad - “ aeT w 


(80) 


For charred surfaces, values of e from 0.7S to 1.0 have been used. 


Homogeneous Combustion Energy Input Rate 

The model that we use for homogeneous combustion contains the assumption 
that this oxidation occurs at a flame front located at the outer edge of the 
boundary layer. The energy released by the reaction increases the enthalpy 
potential and this, in turn, increases the convective heat transfer to the sur- 
face. Thus, only a fraction of the energy released actually goes into the 
body. In the calculation, we make use of the mass rate of air flow per unit 
area into the boundary layer, m a ^ r = (P a ir v air^ 5 as calculated in equation 
(25). We also use the rate of expulsion of pyrolysis gases, = in (0, t) , 
as evaluated in equation (20b) or (21c) for y = 0. 

The energy liberated by combustion of a unit mass of pyrolysis gas, h c , 
can be described as a function of temperature and extent of oxidation which 
will depend on a number of factors (evaluated in the material properties sub- 
program) , or it can be conveniently approximated as a read-in constant. (The 
latter is within the framework of the simple flame front model assumed.) A 
value of about 11,000 cal/g has been used for h g C . If there is a deficiency 
of oxygen, the value of hg C must be corrected to an effective value, h^e, 
which is the actual energy released per unit mass of pyrolysis gas. We 6ave 


l gce 

- h gc 


when 

™ox > 
ex 

0 

(81a) 



/ ™ox \ 

when 

™ox 

ex 

0 

(81b) 

'gee 

- h gC | 

l">ox ) 
\ theo/ 


where m ox and m ox are evaluated in equations (26) and (27), respectively, 
theo 

and m ox is evaluated from equations (28) and (29b) . The rate of energy 
ex 

liberation at the flame front is 


E ff 




(82) 


The increase of enthalpy potential at the edge of the boundary layer 
increment of energy, Eff, per unit mass of external gas entering the 
layer. In km 2 sec -2 , we have 


V 


0.00836E ££ 


m. 


air 


is the 
boundary 


(83) 
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where Aair is evaluated in equation (25) . The heat received by the body 
from the homogeneous combustion is in proportion to the increment in enthalpy 
potential, and this is: 


q cg 


V 2 


q W^\f2 

0.00836CpT w 


(84) 


Heterogeneous Combustion and Erosion Energy Rates 

We assume that all energy liberated by heterogeneous (solid) combustion 
and erosion (material removal by other than combustion or sublimation) is 
received by the body since these are surface reactions. We include both 
energies in a solid combustion term. 


Tcs ^chcomb^scomb + ™chero^sero 


(85) 


The energies, h S Q 0in b and h serQ , can be made functions of temperature (in the 
material properties subprogram); or they can be assigned constant values. The 
combustion energy, h scom ^, is per unit mass of char consumed, and the fraction 
of noncombustible material in the char must be considered in selecting ^s com b 
(see eq. (43)). (If the noncombustible material melts off with a heat of 
fusion, h scom b should be modified to account for this energy.) Values of 
h SCO mb between 1200 and 2500 cal /g have been used. The erosion energy, 
h se ro> is normally assigned the value, zero. 


Sublimation Energy Rate 

The sublimation of char is treated as a surface reaction which removes 
heat from the body. The heat of sublimation of the char, h su ^ , can be set up 
as a function of temperature, but it is most conveniently read in as a con- 
stant positive number. It is an energy per unit mass of char consumed and 
will therefore depend on the composition of the char. One will normally use a 
heat of sublimation for carbon reduced by dividing by Cj^b to account for 
the fraction of carbon in the char. If the other material in the char is con- 
sidered to melt with a heat of fusion when the carbon sublimes, a correction 
for this should be included in h su b . A value of h su b = 2700 cal /g has been 
used for some chars. The energy term due to sublimation turns out to be very 
small unless the flight velocity is extremely large, as from a planetary 
return . Thus 


^sub ^chsub^sub (^) 

where q su b is positive for energy put into the body. 
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Boundary Conditions 

A boundary condition for equation (1) is needed at the front surface of 
the body (y = 0) . For the semi-infinite formulation, this is the only bound- 
ary condition. For the finite depth formulation, a back surface boundary 
condition is also required. 

The front-face boundary condition is obtained from a surface energy 
balance. 


% + ^cg + Qcs + ^sub + ^rad ^sw 


/8t\ 
l 3y/ 


Equation (87a) describes a relationship between (3T/9y) w and the energy rate 
terms on the left side of the equation. The evaluations used for the energy 
terms are listed in preceding sections. The values of the various energy 
terms depend on external flow conditions and on surface temperature; some of 
the evaluations depend on conditions in the interior of the body (including 
past history). The stability of equation (87a) is discussed in the section, 
Some Features of the Computer Program. 

In the finite depth formulation, the heat shield is supported by a back- 
ing material (which may be a composite) . The back surface energy balance for 
the heat-shield material is: 


/ 3T\ 

V y h 


(87b) 


where qgp > 0 means that heat is transferred into the heat shield. For the 
backing material, we also have an equation similar to (87b) . The two surface 
temperature gradients depend on the back- face temperature and on internal tem- 
peratures in the heat shield and the backing material, and therefore they 
depend on the time history of the problem up to time, t. Complicated cases 
can require iteration to obtain n nearly exact" solutions for the back-face 
temperature. In actual numerical calculations, iteration is almost never 
required for good approximations . One simple method is to use K sBF from the 
previous (known) time line. A polynomial passed through the last several 
points of the current (being calculated) time line gives a linear relation 
between q B p and T^p. The same procedure can be used for the front of the 
backing material if it is solved by finite differences. Thus, two relations 
between q^p and T^p provide a solution. 

The back-face temperature, Tgp, and the gradient in the heat shield, 
(3T/3y)gp, are solved for in a back surface subprogram, QBACK. One particu- 
larly simple case (programmed in a QBACK) considers the backing material as a 
heat sink of finite heat capacity at a uniform, time-varying temperature, Tgp. 
For this case 




(87c) 



where c" is the heat capacity per unit area of the backing material. This 
case is also solved by using a polynomial to relate the unknowns, q B p and Tgp, 
in equation (87b). The second relationship, in equation (87c), is obtained 
by approximating (BT/St)^^ as a finite forward difference quotient using the 

unknown, Tgp. The adiabatic case is obtained with c* = 0. 


Trajectory Equations 

For flight cases with aerodynamical ly calculated heating rates we need 
trajectory equations which are solved simultaneously with the other calcula- 
tions. Trajectory equations are not used for flight cases with arbitrary 
heating rates (programmed as inputs) . 

The trajectory equations are in a trajectory subprogram. We use two- 
dimensional trajectory equations for entry in a meridional plane and allow the 
body mass to vary. Equations (88) through (95) are identical in form with 
equations (52) through (59) in reference 1. (See also ref. 9.) 



dD 

dt 

We use a hypersonic approximation that 
constant . 


2 + y 2 

(88a) 


(88b) 

•W •] ' V 

(89a) 


(89b) 

Dv 

' s h 

(90) 


iders the ambient enthalpy as a 


v =V V i + E 12 


(91a) 


E 12 = 0. 008361^ (91b) 

where h^, is the average ambient enthalpy in cal/g and E 12 has the units 
km 2 sec" 2 . For earth entries, 0.5 has been used for E 12 . 
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For earth entries, we generally use an option in the computing program 
that puts in the ARDC atmosphere. This is actually approximated as an exponen- 
tial atmosphere with three successive values of the scale height, S^, which 
automatically change with altitude as programmed. The option also sets 
Rpl = 6440 km and gpj = 980 cm/sec 2 . Another option, for arbitrary planet 

entries, accepts read-in values for Rpj and gpj, arbitrary scale heights 

(initial, two intermediate, and a final value), and atmospheric densities at 
which the scale heights change. 


The quantity, M/CpA (g/cm 2 ) , is needed in the trajectory equations. The 
variation of M/A has been related empirically to the surface recession, X, 
as in reference 1. This is equivalent to assuming a geometry for the 
recession shape. 

x = (x).^ + E ia x + c 3 ( eC4X - *)] C 92 ) 


The variation of the drag coefficient, C D , through the free-molecule, transi- 
tional, and continuum regimes is represented by 


where 


1 

- 


C D = C DC |1 * E, e-lSCRD) C1*E 13 ) 


) 


Eq = 


C DFM " C DC 


•'DC 


(93) 

(94) 


and E 9 depends on body shape. The free-molecule drag coefficient, Cp^, may 

be given the value, 2. The parameter, E 13 , has some dependence on body shape 
(and flight conditions), but will often be assigned the value zero, which is 
the value for a sphere in air. Combining equations (92) and (93) yields 


M 

c d a 



+ E 1 gX + Cg(e E ' 4X - l)J 
1 . E, e -”CRD)CUE 13 )l 


(95) 


In equation (95), C D q is grouped with the initial value of M/Cp^A; this ini- 
tial grouped quantity is read into the computing program. The empiricism in 
equation (95) can also be used to account for any change in Cpc with change 
of body shape. 
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Calculated Quantities of Interest 

The two quantities listed below are calculated by the computing program 
because they are considered to be of interest. They are not needed in any 
other calculations. They are also listed in reference 1. 

Thermal thickness -We define thermal thickness as 


fT - T ) 
1 w o J 

(3T/9y) w 


(96) 


This quantity is meaningful only when the wall temperature gradient is 
negative. 

Stored energy comparison with exponential temperature prof ile -We compare 
each temperature profile with an exponential profile having the same (3T/8y) w . 
We define the quantity, <j>, as the ratio of energy stored to the energy stored 
by the corresponding exponential profile (with form similar to equation (122)) 
Both energies are calculated as constant property approximations. 


! Y BF (T - T 0 )d yi 

♦ = ° (97) 

A ( T w - V 

This quantity, <J>, is only meaningful when A > 0, or the wall temperature 
gradient, (3T/9y) w , is negative. 


Energy Balance 

As an indication of the accuracy of solutions obtained by the numerical 
finite difference method used, an energy rate balance is made from the solu- 
tion at hand. The group of energy rate terms listed below are evaluated and 
summed up; the residual of the sum is the error in the energy rate balance. 
The magnitudes of the energy rates are also of interest since these show the 
dispositions of the various energies. 

The sum of the energy rate terms includes q w , q rac j, q C g, ^cs* anc ^ ^sub’ 

These are evaluated in equations (66), (70), or (73d) for q w , (80), (84), 
(85), and (86), respectively . These are the terms on the left side of the 
front-face boundary condition equation (87a) . The energy rate into the back 
surface, qgp> is given in equation (87b). This term is zero for the 
semi-infinite formulation. 

We require an evaluation of the stored energy in the material. For the 
enthalpy of the solid material (cal /g) , we have, from equation (13b): 
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h S = fZ c 5 dT l = /o [(1 - J)C Ch + JC p] dT l 


(98) 


For the enthalpy of the gas (cal/g), we write 

hg = Cpg dT^ - h sl >ef (99) 

where -h §re f is the energy of formation of the gas at 0°K. (For a pyroly- 
sis reaction that is endothermic at 0°K, h sre f < 0.) The enthalpy per unit 
area per unit depth is then 

h(y, t) = p s h s + Pghg (100) 

where p s and pg are evaluated in 
stored energy is then 

h^ft) = 

For the rate of increase of stored 

dhrp 

9stor = dt 

This is evaluated by the computer program as a finite difference quotient, the 
numerator being the difference between values of Ep for adjacent time lines 
and the denominator being At. With the finite depth formulation, ygp 
shrinks with time. With the semi-infinite formulation, ygp remains constant, 
but actually an infinite depth is being represented by a finite one. This 
finite depth must be sufficiently deep that what is cut off is of no conse- 
quence. This requires that, to within a desired degree of accuracy: 



Ps (yBF» ~ D p 
m g (y BF , t) = 0 


equations (11) and (18) , respectively. The 


y BF 

/ h(yi, t)dyx 
o 


( 101 ) 


energy, 


d r 7BF Cr ^ 
^ / Myi* t ) d yi 


( 102 ) 
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These requirements also must hold in the immediate vicinity of the (ficti- 
tious) back face. The above requirements do not hold for the finite depth 
case . 


The convection energy rate of solid and gaseous ablative material (posi- 
tive out), for the semi-infinite formulation, is: 


Ivcon = v S [ p s^BF’ ^s^BF’ U “ P s (°» ^ + t)h g (0 > ^ 

(103a) 

The depth requirements noted above for a semi-infinite representation insure 
that q vcon is correctly evaluated for this case. 

With the finite depth formulation, the control volume is shrinking and no 
material is entering it; thus: 


^vcon 


= -v s P s (0, t)h s (0, t) 


+ m g (0, t)h g (0, t) 


(103b) 


The energy rate balance is: 

^res + ^rad + ^cg + ^cs + ^sub + ^BF " ^vcon “ ^stor 

where qgp = 0 for the semi-infinite formulation and where q res is the 
residual term in the energy rate balance. We define (arbitrarily) the error 
in the energy rate balance as 


qres 
e — 

when 

lq w l > 1 

(105a) 

1 q w 1 



e = qres 

when 

|q w l * 1 

(105b) 


We also calculate a cumulative energy balance. This balance shows the 
total size of the various energy terms and the error accumulation. We compute 
the integrals: 


Qw = T qw dt l 

H 

Qrad = ^rad d ^l 
Qcg = / Reg 

t A 


(106a) 

(106b) 


(106c) 
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(106d) 


Qcs / ^cs dt l 
1 

Qsub = / Qsub d ^l (106e) 

t i 

QBF = / t QBF dt l (106£) 

*i 

(where Qgp = 0 for the semi-infinite formulation) . 


Qvcon J* ‘Ivcon dt l (106g) 

t i 

Qstor = ^fCt) ” ^rC^i) (106h) 

The accumulated residual is 

Qres = Qw + Qrad + Qcg + Qcs + Qsub + Qbf ~ Qvcon “ Qstor (107) 


Representations o£ Physical Properties 

In the analysis, the quantities, pp and p c ^, the densities of virgin 

plastic and pure char, respectively, are taken as constant. Other properties 
of the virgin plastic and pure char, such as specific heat, enthalpy, and 
thermal conductivity are generally functions of temperature. The temperature 
variation applies also to the specific heat and enthalpy of the pyrolysis gas. 
The temperature-dependent properties are evaluated for each given material in 
a material properties subprogram for the material. The properties data, 
which are in the subprogram, may be in tabular or equation form. 

Some physical properties can be described for a number of materials by 
more or less universally accepted equation forms (e.g., kinetic reaction 
rates), as given in the analysis above; the constants in the equations are 
"read-in 11 quantities. These equations are "built into" the analysis, but more 
complex representations can generally be obtained by prescribing a variation 
of the constants through the use of the material properties subprogram. 

For a given material, the physical properties of the virgin plastic and 
the pure char are generally determined by standard measurement techniques. 

The thermal conductivity of the pure char is probably the most uncertain quan- 
tity. (Some measurements even indicate a hysteresis effect for this quantity; 
see ref. 10.) For a given material, comparisons should be made between abla- 
tion calculations and wind-tunnel measurements, because the "effective" values 
(in the wind tunnel) of the physical properties, particularly K^, may be 
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somewhat different from those values determined from direct property measure- 
ments. When this occurs, the K c ^ is usually adjusted to obtain good wind- 
tunnel comparisons, and the adjusted values are also used for flight 
calculations. In the material properties subprograms, adjusted values of 
properties are usually put in. However, the adjustment can often be made by 
applying a constant multiplier to the measured values. Then, one would enter 
measured values and use an available (open) read- in constant as a multiplier. 


NUMERICAL ANALYSIS AND COMPUTATION PROCEDURES 


In the above analysis, the energy equation (1) with its boundary condi- 
tions, equation (87), is the principal equation system to be solved. Most of 
the other equations listed are used to evaluate quantities that appear in 
equation (1) and its boundary conditions, equation (87). Cases of interest 
are not generally amenable to analytical solutions, so numerical methods must 
be used. Equation (1), a parabolic partial differential equation, is solved 
numerically by finite differences using an explicit (forward difference) 
scheme. As the solution advances from one time line to the next, the first 
quantity solved for is the finite-differenced temperature array, T(M),from 
which the various temperature-dependent physical properties are determined. 
Also, the quantities, K(M) , then z(M), leading to J (M) and p s (M) , are calcu- 
lated by means of the temperature array. At this point, we also calculate the 
rate of generation and flow of pyrolysis gases, w (M) and m (M) , and the sur- 

& o 

face recession rate, X. These quantities become new inputs to equation (1) in 
advancing to the next time line. 


Finite Differencing 


In the finite differencing scheme used to solve equation (1), the partial 
derivatives of the temperature, T, are represented by the following difference 
quotients : 


|- -Ay 


° r 
1 

o J 


(«) 

, _ t M,N+1 - t M,N 

M,N At 

(108) 

(nS 

i _ T M+1 ,N ‘ T M-1 ,n 

(109) 

w 

M,N 2Ay 


9 2 T\ _ t M+1,N 


( 9 ) 

X 7 M,N 


2T m,n + T M-1,N 


(Ay) ' 


( 110 ) 


M - 1 


M 

y 


M +- 1 


where M-l, M, M+l are grid point num- 
bers on the depth (y) scale, and N, 
N+l are numbers on the time (t) scale 
Sketch (e).- Finite differencing as shown in sketch (e) . Finite incre- 
scale. ments are indicated by the A symbol. 
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Stability and Accuracy of the Finite Difference Equation 

In solving a parabolic partial differential equation by forward differ- 
ences, a stability requirement is necessary so that round-off and other errors 
do not grow excessively large. The forward finite difference form of equa- 
tion (1) has an approximate stability parameter requirement 


Z x = 


AtKs 

(Ay) 2 [p s c s + p g c pg ] 


1 

5 2 


( 111 ) 


that may not be sufficient to guarantee stability of the entire calculation. 
If an instability develops (say, at the front face, or at the back face for 
finite depth calculations) , as a general rule one can arbitrarily reduce Z x 
until stability is obtained. This is discussed further in the next section. 
The stability parameter, Z x , is printed out for each finite-differenced grid 
point at each time printed so that the stability of the calculations can be 
monitored continuously. Usually Z x will be largest with the smallest Ay 
used in finite differencing. This is Ay 2 , which is shown in sketch (h) in 
appendix C. 

A gross check on the accuracy of numerical solutions is provided by the 
running energy balance and the cumulative energy balance described above. An 
additional check, standard in numerical work, can be made by adjusting the 
time and space increments, At and Ay, and noting the effect on the numerical 
solutions. This check determines the adequacy of representing the partial 
derivatives by difference quotients with the finite increments as selected. 


Some Features of the Computer Program 

Damping -Within the requirement of the stability parameter, Z x < 1/2, the 
overall calculation may still not be stable mainly because Arrhenius type 
functions of temperature are used that are quite volatile at temperatures 
above a threshold value. In general, the calculations become more stable when 
smaller Z x values are used. This can be accomplished by making At smaller 
or Ay larger. However, if Ay is made too large, the space derivatives are 
poorly represented by the difference quotients used, particularly at the front 
face where the space derivative may be quite large. A value of Ay too large 
may then give an excessively high surface temperature or require an exces- 
sively low temperature at the adjacent finite-differenced grid point. Thus, 
too large a value of Ay can render the front- face boundary condition 
(eq. (87a)) unstable. Reducing the At increment to very small values will 
improve stability but will increase computing machine time, which varies 
inversely as At. 

To avoid excessively small time increments, an option is provided in the 
computing program that allows damping of the calculated changes in certain 
quantities. These are calculated quantities that tend to n overshoot n or oscil- 
late, and damping serves to suppress the oscillation. The damping formula 
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used is: 


X c = X _ i + dfXp - X_ ! ) (112) 

where X is the quantity considered (modified by subscripts); X _ i represents 
the corrected value of X from the previous time line, Xp is the predicted 

value of X for the current time line as obtained from the equations given in 
the analysis, and X c is the corrected value of the quantity, X. The symbol, 
d, represents the damping constant where 0 < d < 1. It is seen that d < 1 
proportionately reduces the change in the quantity, X, from one time line to 
the next. When d = 1, the predicted and corrected values of X are identi- 
cal, and there is no damping. When damping is used, it overrides or modifies 
the normal computations. Damping is available in the computing program to 
correct any or all of the following quantities separately • (y 9 t ) , 

™chcomb> *chero’ ™chsub’ The corrected values are calculated as the step 
immediately arter predicted values are determined so that subsequent calcula- 
tions (and determinations of other quantities) use the corrected values. 
Separate damping constants are read into the computing program for each of the 
quantities listed; each damping constant can be changed twice during a run. 

The latter provision allows a discrete application of damping when it is 
needed, as for example, during rapid changes in external conditions acting on 
an ablator. 

The application of damping should be minimized (d 1) to that amount 
necessary to smooth or reduce oscillations. Excessive damping can delay the 
variation of a calculated quantity in its approach to the correct calculated 
value. The use of damping requires some judgment, and may involve some trial 
and error and repeated calculations. In this connection, important criteria 
are the residual rate, q res > and the accumulated residual, Q res , in the energy 
balances. The residuals should, obviously, be minimized. In damping the T w 
calculation, another important criterion is the change in (3T/3y) w ; the method 
of using this criterion is explained below. 

Front-face bound ary condition -The front-face boundary condition equa- 
tion (87a) furnishes a nonlinear relationship between T w and (3T/3 y) w , but 
also involved in this equation are the internal temperature distribution and 
the extent of degradation of the material. The terms q w , q C g, and q cs are 
affected by the rate of expulsion of pyrolysis gases, and the rate of gas 
expulsion depends on the interior temperature distribution and the extent of 
degradation. The extent of degradation at a time, t, obviously depends on the 
entire past history of the problem. The energy rate terms on the left side of 
equation (87a) also depend on external conditions which may be changing with 
time (as for example, with an entry flight). We see that, in order to solve 
equation (87a) for the boundary conditions at time, t, we must have solved the 
entire problem up to time t. 

Using a forward difference scheme to solve the partial differential equa- 
tion (1), we obtain the internal temperatures for the next time line. For the 
front surface temperature, we seek a numerical solution to equation (87a) . We 
could iterate a solution to equation (87a) by requiring that T w and (3T/3y) w 
fit a polynomial that passes through several interior grid points. This 


33 



consumes computing time, because terms on the left side of equation (87a) 
depend on T w , and would have to be recalculated with every iteration cycle. 

In place of iteration, we use (with little loss of accuracy) the follow- 
ing scheme in the computing program. We set up a two-term Taylor series: 



(113) 


where the subscript -1 means from the previous time line, and 




Equation (114) contains (d/dT w )q,I s . The term q cs as evaluated in equa- 
tion (85), includes the heat liberated by erosion (if h sero i 0). For equa- 
tion (114), q cs is modified and q£ s is defined as: 


^cs ^chcomb^scomb 

thus the erosion term is neglected in calculating the predicted surface tem- 
perature, T W p. This is exact in the usual case when h sero = 0. Otherwise, 
it is approximate, but an erosion term in equation (114) should always be 
small. To include this term would involve several equation forms to corre- 
spond to the different empirical forms for evaluating erosion itself 
(eq. (46)). Equation (113) has the form: 



A + B(T W p - T w _ j ) 


(116) 


We also evaluate (9T/9y) by means of a cubic passed through the first four 

wc ' 

(finite difference) calculated grid points of T. (See sketch (h) , 
appendix C.) 
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The increment, Ay 2 , is the forwardmost Ay used in the finite difference 
equation, as indicated in sketch (h) . Equation (117) has the form: 



11 

6Ay 2 



+ C 


(118) 


Combined equations (116) and (118) become: 


C - A * BT W.! 
T = i 

Wp B + (11/ 6Ay 2 D 


(119) 


The computer prints T^ as a quantity of interest. It is inserted into 
equation (112) to obtain the corrected surface temperature, ^wcj which is the 
value used (and also printed out) . We obtain a corrected calculated value of 
the front surface temperature derivative by rewriting equation (117) . 




jLL r + 3T 
6 [ wc + ^l2 


— T + 
2 L3 


It 1 

3 L4j 


( 120 ) 


Once a corrected value of surface temperature, T wc , is determined for the 
new time line, the terms on the left side of equation (87a) can be evaluated. 
Then, from equation (87a), the surface temperature gradient, now labeled 
(8T/8y) we , that satisfies the surface energy balance can be evaluated. If the 
calculations would be perfect, we would have 

■ a. 

or the two methods of calculation should yield the same answer. (Iteration 
would accomplish this to a desired degree of accuracy.) For comparison, the 
computer prints both evaluations of (3T/3y) w . The closeness of the comparison 
furnishes a criterion for possible adjustment of the time and distance spacing 
and, if necessary, damping in the T w calculation. 

The evaluation (8T /3y) we is used as (3T/3y) we in equation (113) in 

making the calculations for the next time line. With the procedure used, a 
small oscillation is expected in the value of T w about the true value, but 
the oscillation is self-correcting. For example, if T w is too high, the 
energy rate terms on the left side of equation (87a) will be too small, and 
the absolute value of (3T/3y) we calculated from equation (87a) will also be 
too small. This will reduce the value of Tw as calculated for the next time 
line, so that corrective action will occur from one time line to the next. 

Back boundary condition- In the semi- infinite formulation of the problem, 
there is no back boundary condition. In the numerical solution of 


35 


equation (1), however, we are actually representing an infinite depth with a 
finite depth having a finite number of grid points. The deepest point is 
assigned the temperature, T 0 , which it holds for all time. An initial tempera- 
ture profile must be selected such that the grid point adjacent to the deepest 
point has the temperature T 0 to within a desired accuracy. Also there 
should be no degradation at either of these points (J = 1). The point adja- 
cent to the deepest point is watched during a calculated run. It should expe- 
rience essentially no degradation (J = 1) and the temperature should remain at 
T 0 to within some desired accuracy. If these criteria are not met, more 
depth is needed (to represent infinity) . 


In the finite depth problem, we have the back boundary condition, equa- 
tion (87b) . This is solved in a back surface subprogram for Tgp and 
(3T/9y)gp, and qgp is also obtained (in the main program) from equa- 
tion (87b) . For tne simple finite heat sink type of backing, the calculation 
also makes use of equation (87c) . 


As the front surface recedes 
with time, the material depth 
decreases. Since the origin of the 
coordinate system is kept at the 
front face, the spacing between the 
back-face grid point and the adja- 
cent point is allowed to decrease. 
(The printout shows the depth back 
of the front face of each grid 
point.) As indicated in sketch (f ) , 
the initial array consists of MF 
grid points with an MF+1 point that 
coincides with the MF point at the 
back face. As the back face of the 
material moves forward, the MF+1 
point moves with it. The MF point 
remains fixed in space, and thus becomes a virtual point. The temperature of 
the virtual point, T^p, is calculated as a linear extension of the profile in 
the material. This temperature is used in calculating the next time line 
(because of the convenience of using the uniform spacing of the virtual 
point). It is printed out as a quantity of interest, not as part of the regu- 
lar array of y, T, J, etc., at the grid points. When the surface recession, 
X, exceeds the spacing between grid points at the back (in the sketch the back 
face, BF, moves ahead of grid point, MF-1), one n'ode is dropped at the back. 
Thus, the numerical values of all the grid points shown in the sketch are 
reduced by one, and the calculation procedure is the same as before. 


MF-2 


MF- 1 


MF + I 
MF 


BF 


MF-2 MF-I MF + I MF 

% * (•) 


BF 


at time, t 


Sketch (£).- Back-face grid spacing. 


At the back face, it is possible for instability to develop when the sta- 
bility parameter, Z x < 0.5. If this occurs, Z x should be made arbitrarily 
smaller until stability is obtained. 
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Starting Values 


To start the solution to equation (1) (which is parabolic) , it is neces- 
sary to assign initial conditions. The particular selection of initial condi- 
tions is generally not critical as their influence damps out in a short time. 
The equations developed in this section are similar to those developed in 
appendix B of reference 1 . 

Continuatio n type program- With a continuation type computing program, 
none of the equations shown below in this section are used. The initial tem- 
perature profile is read into the program (also the profiles of extent of deg- 
radation and rate of gas expulsion) . These arbitrarily read-in profiles are 
generally the final profiles of a calculation that was interrupted and is now 
being continued. 

Starti ng type programs- With starting type computing programs, we assume 
we have virgin plastic material (J = 1 everywhere) . The temperatures should 
be sufficiently low that they could exist before the onset of ablation. The 
initial temperature profile that we assume is exponential. 


y 



T i (y) = To + CT wi - T 0 )e Ai 

(122) 

If T w i is selected near T 0 , this profile is a small perturbation on a con- 
stant T 0 profile. We take the y derivative of T^ at y = 0, equate it to 
the ratio of the heat flux (eq. (87a)) to thermal conductivity, and solve for 
the initial thermal thickness, A^. (See eq. (96).) 


^pwi ( T wi ~ T o) 

i ^wri 

(123) 

where 

Hwri "" ^wi + ^radi 

(124) 


is the combined initial convective and radiative heat flux. 

Arbitrary heating-rate programs (s tarting type) -Equations (124), (123), 
and (122) are used in that order to obtain the starting profile. The quanti- 
ties, q w ^ and q ra di> in equation (124) are obtained from programmed inputs. 
The initial convection rate, q w i > is the same as q G obtained from an inver- 
sion of equation (72), where q ocw is a programmed input and V 2 (eq. (56a)) 
is obtained from h st , which is also a programmed input. For the initial net 
radiation input, q ra ^i> re-radiation is neglected (eq. (80)), and q ra di is 
equated to q R , which is a programmed input. 

Calculated heating-rate programs (sta rting type) - Equations (124), (123), 
and (122) are used in that order for the starting profile. However, as 
described under Flight Cases, one option allows Ai to be assigned a value, 
instead of obtaining it from equation (123) . The initial convective and 
radiative heat fluxes needed in equation (124) are obtained differently for 
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wind-tunnel and flight cases, as described later. Calculating q w ^ and q ra ^^ 
requires the initial free-stream density, D^. For wind-tunnel calculations, 
will be known; for flight, Di can be assigned or calculated as will be 
described. 

Ablation in a w ind tunnel-The starting of a wind tunnel is considered as 
a sudden step of a heat flux. We assume that i|/ = 1, so 

q wi = %ci 025) 


and 


^wri ^oci + ^radi (126) 

We evaluate q oc i from equation (62) and q ra di from equations (79) and 
(80), assuming negligible re-radiation and using T w ^ for T w . If T w ^ is 
small, q wr i is approximately constant for a short period of time (eq. (62)). 
We approximate this case as the classical conduction problem with constant 
heat flux and constant properties (refs. 1 and 11), and we calculate the time 
at which the front face arrives at the assumed temperature, T w ^ . 

C T wi - V^PpViV'i 

t- = “ — (127) 

4q 2 . 
n wri 

We can set T w ^ > T 0 ; the greater value is not necessary, but it gives the 
computing program a smooth start. 

Flight with arbitrary initial conditions-In starting flight calculations, 
we can use an assigned initial velocity, V^, and a flight-path angle, y^, at 

(arbitrary) time, t^ = 0. Equivalent to a starting altitude, the initial 
atmospheric density, D-[, can be assigned as well as the initial thermal thick- 
ness, A value for T w ^ is assumed and the initial temperature profile is 

obtained from equation (122) . Using this procedure we do not require that Ai 
be consistent with the relationship in equation (123) , but we still calculate 

°lwi 5 ^radi * anc ^ ^wri * 

Initial condit ions for entry flight-For entry flights an alternative 
starting procedure is described in appendix B of reference 1. An assigned 

and y^ and an assumed T w ^ are used to calculate the entry into an expo- 
nential atmosphere. Heating during this entry raises the front-face tempera- 
ture from T 0 to the assumed Twi- The convective heating rate during this 
initial entry phase is considered to be a free-molecule type so we have: 


and Kt u 


q wi ~ q FMi K tu 0 2 8) 

in equation (61) specializes to K. Thus 

q wi = Kq FMi C 129 ) 
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For the initial combined convective and radiative heating rate, equation (124) 
is written as 


^wri ^FMi + ^radi 


(130) 


We evaluate f rom equation (64) and q ra dj_ from equations (79) and (80) , 

and neglecting re-radiation, we have the heating rate at time, t^ = 0. 


KA^D.V; 


Ec E, 


<W - o?0836 1 (V i ‘ <>. OOSSOT^j) * KE^D^y!: 6 


(131) 


where we use as an approximation for V*,. . 

atmosphere, we write the density as 

D = c e 



In the assumed exponential 


(132) 


where Alt represents altitude. During the entry phase from t = to 
t = tj_ = 0, the enthalpy potential can be approximated as constant, and the 
heating rate, q wr , will increase exponentially as D increases. We can 
rewrite equation (131) for q wr and substitute the variable D for , and we 
can integrate the equation from t = -« to t = tj = 0 for the total heat 
absorbed. 


Qtotal 


Sh.KA^DiCV? - 0 . 00836c p T wi ) 

+ 

0.0836 sin 


Sh.KE 4 R.D? 5 V? 6_1 
1 _ l i i 

E 5 | sin y.| 


(133) 


The total heat absorbed is approximated as 


JMF 

Qtotal = Pp c pwi ^ (T " T ° )dy 

and using the profile equation (122) we have approximately 


(134) 


Qtotal * Pp c pwi( T wi “ T o) A i 


(135) 


We substitute q wr i from equation (131) into equation (123) for A^, and put 
A-l into equation (135) to obtain: 


Qtotal 


Pp c pwi K pwi ( T wi 


To) : 
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(136) 
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We eliminate Q-^otal ^y combining equations (133) and (136) and we obtain an 
equation of the form 



(137) 


where K a and are constants and is the only unknown. 


The procedure in starting an entry calculation is to assume a T w -[, and, 
knowing the scale height far out in the atmosphere, S}^, calculate from 

equation (137) . Next obtain q wr i from equation (131) , and from equa- 

tion (123); then put A^ in the profile equation (122). The assumption of 
the T w ^ determines the Di, or equivalently, the altitude at which to start 
time zero. In order to have a finite starting altitude, it is necessary to 
select T wi > T Q . 


Use of Programs in Various Approximations 

The analysis outlined above is for a blunt-body axisymmetric stagnation 
point. The approximate use of the analysis in the stagnation region of an 
arbitrary blunt body is considered quite valid. The analysis is actually one- 
dimensional, but in the stagnation region, lateral gradients of such quanti- 
ties as temperature (and resulting heat transfer) are very small relative to 
gradients normal to the surface. The heat loading for any point for which 
calculations are being made must be known or calculated. The use of the anal- 
ysis for points far fron the stagnation region is still more approximate, but 
it is often acceptable when these points have relatively small heat loadings, 
as is often the case. 


In the analysis, the stagnation pressure, p t or p s , appears a number of 

times in the development. This should be considered to be the local external 
pressure, p w , for off-stagnation calculations. This applies to equations (15), 
(24b), (37), (39), (40), (46a), (50), (51), and (52). For a flight case with 
calculated heat loadings, the constant, Aj , in equation (54) is adjusted to 
give the local external pressure p w . For the calculations of local convec- 
tive and radiative heat loadings, the constants to be adjusted are A 4 , A CQ , 

A 5 , and E 4 in equations (62), (64), (77a), and (79), respectively. In equa- 
tion (62), if desired, the effective nose radius, R, may be changed along with 
A 4 by adjusting the constants in equation (60) . For calculation of the oxygen 
available for combustion off the stagnation point, the constant in equa- 
tion (24b) may be assigned values other than unity. These same adjustments 
are also made for wind-tunnel cases with calculated loadings, except that the 
local pressure, p w , is read in (FORTRAN symbol PT2) . For off-stagnation calcu- 
lations with arbitrary heating rates, the constants to be adjusted are <jT, A*, 
and A 5 in equations (24b), (54), and (77a), respectively. For these cases, 
equation (54) is inverted to yield the free-stream density, D, from the pro- 
grammed p t or p w . The radius, R, is programmed as a function of time. 
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although R may be a fictitious effective quantity for many of these cases. 
The program calculates and uses (but does not print) A 4 from an inversion of 
equation (62) . 

The stagnation-point analysis developed implies the assumption of a lami- 
nar boundary layer with a laminar convective heat-transfer rate and diffusion 
rate. An off-stagnation point may have a turbulent boundary layer that will 
introduce a further approximation (and complication) in the analysis and com- 
puting programs. Equation (62) for q QC is considered to be a laminar stag- 
nation point form. A calculated turbulent q QC rate may be obtained by 
adjustment of A 4 (and possibly R) in equation (62) . This rate should be 
satisfactory for the constant conditions of a wind-tunnel calculation. For a 
flight case with calculated heating rates, the evaluation is more difficult. 

It may then be necessary to program A 4 (and possibly R) as functions of V 
and D in the atmosphere considered, to account for effects of time- varying 
Mach and Reynolds numbers for the point considered. (This could be done in 
the material properties subprogram.) For arbitrary heating-rate calculations 
the situation is somewhat better, as q ocw is known and q Q can be obtained 
from equation (72). An effective R, as a function of time, can be programmed 
to be consistent with equation (62) . For both the calculated and the arbi- 
trary heating-rate types of calculations, A 5 in equation (77a) should^ be 
selected to give a realistic turbulent ip function. The value of <p in 
equation (24b) should probably also be adjusted (or programmed) for a 
turbulent boundary layer. 


Illustrative Example 

An example is given to illustrate a typical use of the numerical comput- 
ing programs. This example is typical of Apollo flight conditions. The fig- 
ures given show the time variation of some of the quantities of interest; a 
number of other quantities, not plotted, are also computed and printed out by 
the program. See appendix C for the read outs (for another case) from a 
typical printing subprogram. 


This example calculation is for a body point in the stagnation region 
near but not at the stagnation point. This is then an approximation in the 
use of the program. (Experience has indicated that this approximation is 
quite good.) The trajectory and heat loadings in the example are assigned. 


Figure 1 shows the assigned trajectory while figure 2 gives the total 
enthalpy, h st , and surface pressure, p w . The total enthalpy was calculated by 
combining equation (91a) with the inversion of equation (56) . The surface 
pressure, p w , was obtained by first calculating p t using the hypersonic 


approximation, equation (54), with A^ = 1. 


Then 


p t 2 


was multiplied by a 


factor to obtain p w . 


The multiplying factor 


Pw/Pt 2 


was obtained from the 


NASA Manned Spacecraft Center in a private communication. It had been deter- 
mined from wind-tunnel tests of the model configuration at various angles of 
attack. 
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Figure 1.— Trajectory for illustrative example. 
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Figure 2 — Total enthalpy, h st , and surface pressure, 
p w , at a near-stagnation body point. 


Figure 3 shows the assigned 
cold-wall convective and the radia- 
tive heat loadings as functions of 
time for the body point being cal- 
culated. These loadings obtained 
from the private communication 
noted above had been evaluated by 
using the equivalents of equa- 
tions (62), (72), and (79). The 
heating-rate distributions over 
the surface were mainly based on 
wind-tunnel tests of the model at 
angles of attack. 

The calculated heat-shield 
response to the assigned loadings 
is illustrated in figures 4 
through 7. In figure 4, typical 
calculated results show the varia- 
tion with time of surface reces- 
sion, X , and the sum of surface 
recession and depths of degrada- 
tion of material (to 90 percent 
degraded, J = 0.10, and also fur- 
ther back to only 10-percent 
degraded, J = 0.90). In the lat- 
ter stages of the flight, the heat 
loading has decreased (fig. 3) and 
the surface recession has virtu- 
ally ended, but interior degrada- 
tion of material continues at a 
slow rate because of the slow fall- 
ing off of the internal 
temperature . 



Figure 5 shows the rate of 
mass loss due to pyrolysis, iiig W , 
and that due to surface recession, 

^surface = *chcomb + *chero + *chsub * 
Most ,of the material lost in this 
flight is due to pyrolysis. The 
loss proportion between pyrolysis 
and surface loss is partly deter- 
mined by the properties of the 
material; in particular, it is 
known that with degradation some 
materials produce more char than 
other ablating materials. 


Figure 3.— Cold-wall convective, q ocw , and radiative, Figure 6 shows the calculated 

qR, heat loadings at a near-stagnation body point. values of temperature at the 
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Figure 4.— Calculated heat-shield response; surface 
recession and depth of degradation. 
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Figure 5.— Calculated heat -shield response; 
mass loss rates. 


Typical energy rate dispositions 
at a near- stagnation body point 
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Figure 7.- Calculated energy rate dispositions. 


surface and at several fixed points in the interior of the material. There is 
a time lag of temperature rise with depth in the material , as would be 
expected. During the latter stages of the flight, the surface is actually 
undergoing cooling, and the highest temperatures are then in the interior. 

The curve for a depth of 0.25 cm ends when surface recession reaches that 
depth. For any given time, one can use figure 4 to determine the state of 
degradation at each of the fixed depths for which temperatures are plotted in 
figure 6. 

Figure 7 shows typical energy-rate dispositions for three selected times 
in the flight. At 50 seconds, the loading is moderate (and increasing), as in 
figure 3. Most of the energy loading is accounted for by blockage. Pyrolysis 
is occurring (fig. 5), but the surface temperature is still low (fig. 6), so 
there is relatively little radiation. At 100 seconds, the loading is near the 
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maximum (fig- 3), and most of the energy is accounted for by re-radiation 
because the surface temperature is high (fig. 6) . Blockage accounts for about 
one-third of the energy input, a relative decrease from that at 50 seconds. 

In energy units, cal cm" 2 sec" , however, this blockage is about 50 percent 
higher than that at 50 seconds; the pyrolysis rate is greater (fig. 5), and 
there is more convective heat loading to be blocked. At 500 seconds the rela- 
tive energy dispositions are only slightly shifted from the proportions at 100 
seconds. The main effect here is that the loading has become small (fig. 3). 
The pyrolysis rate is low (fig. 5), so that the blocked energy magnitude is 
small, and the re-radiation is also small because of the low surface 
temperature (fig. 6) . 


An approximate measure of the ablator efficiency is the smallness of the 
absorbed energy term as this is the only energy that actually gets into the 
ablator. This energy is mainly accounted for by the pyrolysis, heating up, 
and loss of gaseous material, and to a lesser extent by the heating up and 
loss of solid material from the surface. A very small amount of the energy 
absorbed is also accounted for by storage of sensible heat in the ablator. 
(Over the entire flight the total sum of this energy becomes almost negligible 
because of the cooling at the end of the flight.) 


A sublimation energy term was (routinely) calculated, but it is negligi- 
ble for this flight because of the range of surface temperatures. Sublimation 
with a charring ablator can be important, and even dominant, for higher speed 
flights such as from planetary returns. Sublimation involves a loss of char- 
ring material at the surface with the removal of energy at the surface. 


DISCUSSION AND CONCLUSIONS 


A general method has been presented for determining the response in the 
stagnation region of the charring ablator type of heat shield. The analysis 
given is actually for an axisymmetric blunt-body stagnation point, but calcu- 
lations obtained from the analysis form a good approximation for the stagna- 
tion region of any arbitrary blunt body when the heat loading can be evaluated 
for the point at which the calculations are made. An even more approximate 
use of the calculations can be made for body points far from the stagnation 
region. This is considered acceptable when these points have relatively low 
loadings. The calculations can be performed for a variety of conditions, such 
as for wind-tunnel tests or flight cases; the heat loadings can be calculated 
by the computing programs, or the loadings can be arbitrarily assigned as func- 
tions of time. Properties of the heat-shield material can be arbitrarily 
assigned. 

The important use of the analysis and the associated computing programs is 
to calculate the performance of heat shields for flight. For any given heat- 
shield material, it is highly desirable to have comparisons of calculated heat- 
shield performance with results from a range of experimental tests in wind 
tunnels. Comparisons with experiment are important, because the effective 
physical and thermodynamic properties of a material that is being used as a 
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charring ablator may not be completely known. When these comparisons are 
good, one can have confidence that the properties of the heat-shield material 
are being adequately represented in the calculations, and predictions of heat- 
shield performance for flight cases should be satisfactory. 

Since heat shields have a finite depth, the most rigorous method of cal- 
culation is the finite depth formulation with an appropriate back boundary 
condition. The semi-infinite type of calculation is faster and simpler. It 
is virtually as accurate if there is a safety factor of sufficient material 
such that there is virgin plastic remaining at the conclusion of the flight or 
test. This is the case because the events of importance occur near the front 
surface where temperatures can be an order of magnitude higher than at the 
back. One can generally say that, if it is necessary to use a finite depth 
calculation to represent the physical situation, the initial depth or initial 
amount of heat-shield material is probably marginal. 

Numerical calculations are carried out by means of a family of computing 
programs. The selection of programs and the options within the programs are 
described in appendix C. This appendix gives detailed instructions in the use 
of the programs, and it contains a sample case for illustration. The programs 
can be obtained from COSMIC, University of Georgia, Athens, Georgia, 30601. 

Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 94035, June 22, 1970 
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APPENDIX A 


PRINCIPAL NOMENCLATURE 


In performing computer calculations, some purely FORTRAN quantities are 
used in the input data which have no counterpart among the symbols listed 
below. These quantities are in appendix C, wherein all FORTRAN quantities are 
listed and identified. 


A 


frontal area of vehicle, cm 2 


A cq> A cv 

Ai 

A 3 

a 4 

a 5 


free-molecule accommodation coefficients for heat transfer and 
sublimation mass loss, respectively, dimensionless 

constant defined by equation (54) 

constant defined by equation (75) 

constant defined by equation (62) 

constant defined by equation (77a) 


B 


blowing parameter defined by equation (76a) , dimensionless 


c ch> c p> c s specific heat of char, virgin plastic, solid material (virgin 
plastic to char), respectively, cal/g °K (see eq. (13)) 


‘rg 


c p 


specific heat of pyrolysis gas at constant pressure, cal/g °K 

heat capacity per unit area of backing material, cal/cm 2 °K 
(eq. (87c)) 

average specific heat, external gas, cal/g °K; defined by 
equation (57a) 


C D^DC^DFM drag coefficient, dimensionless; drag coefficient, continuum; 
drag coefficient, free molecule 


C E 


Arrhenius activation temperature for surface combustion, °K 
(eq. (35)) 


Cj modification factor for thermal conductivity in pyrolysis zone, 

dimensionless (eq. (14)) 

Ckdb reciprocal of mass ratio of carbon/char, dimensionless 

(eq. (43)) 


C ox mass fraction of oxygen in ambient gas, dimensionless (eq. (26)) 
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^oxf 

c oxw 

C x 

3 

Ci,C 2 

C 3 y C 4 

C 6 

d 

D 

e 

e 2 

e 3 

E4.E 5 ,E 6 

Ey 

E 8 
e 9 
E 10 
E 1 2 

E 1 3 

E 1 6 
E 1 7 


mass fraction of oxygen (dimensionless) at outer edge of 
boundary layer after depletion by homogeneous combustion 
(eq. (30)) ' 

mass fraction of oxygen (dimensionless) at the wall (eqs. (31)- 
(33), (36) - (41) ) 

stoichiometric ratio of oxygen to carbon (by weight) in the sur- 
face combustion reaction, dimensionless (eqs. (43) -(45)) 

constants defined in equation (77) 

constants in nose radius (eq. (60)) 

constants in M/A (eqs. (92), (95)) 

constant vorticity correction in equation (62) 

damping coefficient (general), dimensionless (eq. (112)) 

free-stream density, g/m 3 

error in energy rate balance, defined by equation (105) 

stoichiometric ratio (by weight) of oxygen/pyrolysis gas for the 
homogeneous combustion reaction assumed, dimensionless 
(eq. (27)) 

constant in average specific heat of external gas (eq. (57b)) 

constants in gas-cap radiation (eq. (79)) 

constant to account for shift of sublimation equilibrium 
(eq. (50)) 

constant in vapor pressure (eq. (49)) 
constant defined by equation (94) 

constant in average specific heat of external gas (eq. (57b)) 

average ambient enthalpy, km 2 /sec 2 (eq. (91)) 

constant depending on body shape and flight conditions in the 
drag calculation (eqs. (93), (95)) 

constant used in tumbling correction (eq. (61)) 

constant in nose radius (eq. (60)) 
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E 18 
E 1 9 
E ff 

E e » 


§pl 

h 



hgce 

comb 



h sref 


h sub 


constant in M/A (eqs. (92), (95)) 

constant in vapor pressure (eq. (49)) 

rate of energy liberation at flame front from homogeneous combustion, 
cal/cm 2 sec (eq. (82)) 

dimensionless stream function at outer edge of boundary layer and at 
wall, respectively (see appendix B, eqs. (B1)-(B7)) 

gravitational acceleration of planet, cm/sec 2 (eq. (89b)) 

enthalpy, cal/g 

enthalpy of pyrolysis gas, cal/g (eq. (99)) 

energy liberated by combustion of a unit mass of pyrolysis gas, cal/g 

effective energy of combustion per unit mass of pyrolysis gas, cal/g 
(defined by eq. (81)) 

enthalpy of solid material, cal/g (eq. (98)) 

energy liberated by combustion of a unit mass of char, cal/g 
(eq. (85)) 

energy liberated by removal of a unit mass of char by means other than 
combustion or sublimation (called erosion), cal/g (eq. (85)) 

energy released by pyrolysis reaction per unit mass of gas produced, 
cal/g; negative for endothermic reactions; used in equation (1); 
generally a function of temperature, and can be approximated as 
h sr = (Pphp - Pch^ch)/(Pp ~ Pch) “ h g ; can also be approximated as 

a constant; see h sre p below 

a constant; same as h sr if h sr is approximated as constant; also 
used to evaluate hg (eq. (99)); when h sr is made temperature 
dependent as in equation (99) , h 5re f is a reference h sr at a 
reference temperature (usually CrK) ; FORTRAN symbol, HSREF; cal/g 

stagnation enthalpy, cal/g 

energy absorbed by removal of a unit mass of char by sublimation 
(positive value), cal/g (eq. (86)) 


h^ free-stream enthalpy, cal/g 

h enthalpy of interior material, solid and gaseous, per unit area per 

unit depth, cal/cm 3 (eq. (100)) 
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h.'j’ 

j 

J 


K 


K ox 


stored energy (enthalpy) of interior material, solid and 
gaseous, per unit area, cal/cm 2 (eq. (101)) 

order of surface combustion reaction, dimensionless 
(eq. (36)) 

measure of degradation defined by equation (7), dimension- 
less; J = 1 represents virgin plastic; J = 0 represents 
pure char 

Arrhenius rate for pyrolysis reaction; defined in 
equations (4) , (5) 

oxygen diffusion rate constant; defined in equation (32) 

Arrhenius reaction rate of oxygen consumption in surface 
combustion; defined by equation (35) 


K r 

K re 

Ki 

K 2 

K 


K-i 

Kch>V ? s 


Le 


Arrhenius frequency factor (constant) in K ox (eq. (35)) 

defined by equation (40) 

constant, read in (CK1) ; defined by equation (38) (see also 
eqs. (39), (40)) 

constant, read in (CK2) ; modifies char erosion rate, 

dimensionless (see eq. (46) and following discussion) 

fractional time that calculated point is initially exposed 
to approximately stagnation conditions, dimensionless; 
for tumbling correction, equation (61) 

heat-transfer tumbling correction, dimensionless (eq. (61)) 

same as K, but for th_e preceding time line and spaced 
according to the y _ \ array (used in eq. (12)) 

coefficient of thermal conductivity of char, virgin plastic, 
solid material (virgin plastic to char) , respectively, 
cal/cm sec °K (eq. (14)) 

Lewis number, dimensionless 


L 

Dr 

m 


m c 


lift/drag ratio, dimensionless 

exponent in erosion rate calculation, dimensionless 
(eq. (46)) 

rate of mass loss that contributes to convective heat 
blockage, g/cm 2 sec (eq. (75)); (> 0) 
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m chcomb * m chcomb 9 
d max 

m chcomb 
r max 


rate of char loss by combustion, maximum rate by oxygen 
diffusion, maximum chemical reaction rate, respectively; 
g/cm 2 sec (eqs . (43), (44), (45)); (> 0) 


^chero rate of char loss by means other than combustion or sublima- 

tion (called erosion), g/cm 2 sec (eq. (46)); (> 0) 


1 "chsub’™chsFM , rate of char loss by sublimation, free-molecule rate, 

diffusion rate, respectively; g/cm 2 sec (eqs. (53), (48), 
m chsd (52)); (> 0) 


m„ 


m gw 

m surface 

VXj 

M 


M 


cv 




M. 


My 

n 

P 


Ps 


2 


Pt 


2 


Pve 


Pvm 


rate of flow of pyrolysis gas past a point in the material, 
g/cm 2 sec (eqs. (20b), (21c)) (> 0) 

ihg at the surface (pyrolysis gas expelled), g/cm 2 sec 
“chcomb + A chero + A chsub » S /cm2sec 

total mass loss rate of material, g/cm 2 sec (eqs. (74a, b)) 

mass of vehicle, g; also grid point number in finite 
difference computation 

molecular weight of char sublimation vapor, g/mole 
(eqs. (47), (48), (52)) 

molecular weight of pyrolysis gas, g/mole (eq. (15)) 

molecular weight of gaseous products after heterogeneous 
combustion, g/mole (eq. (77a)) 

total mass loss of material, g/cm 2 (eq. (74c)) 

order of pyrolysis reaction, dimensionless (eq. (4)) 

pressure, dynes/cm 2 or atm, as specified 

external pressure, dynes/cm 2 (eq. (IS)) 

pressure downstream of normal shock, atm 

equilibrium vapor pressure of char sublimation vapor, atm 
(eq. (49)) 

modified equilibrium vapor pressure of char sublimation 
vapor, atm (eq. (50)) 


p^ external pressure, atm 
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P ratio of pressure of external gas to modified equilibrium vapor 

pressure of char sublimation vapor, dimensionless (eq. (51)) 

Por porosity of char, dimensionless (eq. (17c)) 

q heat-transfer rate, cal/cm 2 sec; all surface heat-transfer terms are 

positive into the material 

The following are front surface convective heat-transfer rates, cal/cm 2 sec: 

q FM free-molecule rate (eq. (64)) 

q ocw cold wall rate (eq. (72)) 

q oc continuum rate with no blowing (eq. (62)) 

q OQ no blowing rate, bridged between qp^ and q oc (eq. (67)) 

q Q q 0Q corrected for tumbling or oscillation (eq. (68)) 

q^ c continuum rate with blowing (eq. (63)) 

q^ w with blowing, bridged between q FM and q^ c (eqs . (65), (69b)) 

q w q^ w corrected for tumbling or oscillation (eqs. (66), (70)) 


The following are surface heat-transfer rates (other than front surface 
convective) , cal/cm 2 sec: 

q DC . heat-transfer rate into back surface (usually by conduction) 
(eq. (87b)) 


q cg 

q cs 

q cs 

q R 

q rad 

^sub 

9wri 


front surface heat -transfer rate due to homogeneous combustion 
(eq. (84)) 

front surface heat -transfer rate from heterogeneous combustion and 
erosion (eq. (85)) 

front surface heat-transfer rate due to heterogeneous combustion only 
(eq. (115)) 

gas cap radiation rate into front surface (eq. (79)) 

net radiative heating rate into front surface (eq. (80)) 

heating rate into front surface due to sublimation (eq. (86)) 

combined initial convective and radiative heating rate into front 
surface (eq. (124)) 
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The following energy rate terms are not surface heat-transfer rates. They 
appear in the energy rate balance equation (104) , along with the appropriate 
surface heat-transfer rate terms; cal/cm 2 sec. 


^res 


residual in energy rate balance (eq. (104)) 


°lstor 


rate of increase of stored energy in the material (eq. (102)) 


q vcon convective energy rate of solid and gaseous material, 

positive out of control volume; it is the energy taken in 
to heat up and degrade the material (eq. (103)) 


The following energy terms are time integrals of the corresponding q values; 
they appear in the cumulative energy balance equation (107) . 

Qw^BF’^cg* I 

Qcs’Qrad’Qsub’ ! cal/cm 2 (eq. (106)) 

Qstor’Qvcon I 


Qres 

R 

D . 

^air 


V 

S h 

s hi 


s h > S h > S h 
n l u 2 3 


residual in cumulative energy balance, cal/cm 2 (eq. (107)) 

effective nose radius, cm (may be fictitious for off- 
stagnation approximate calculations) 

gas constant for external (ambient) gas, atm cm 3 /g °K 
(used only in eqs . (37), (38)) 

universal gas constant, ergs/mole °K 

planet radius, km (eq. (89)) 

atmosphere scale height, km (eq. (90)) 

initial scale height for entry into an arbitrary atmosphere, 
km (eq. (133) and contained in eq. (137)) 

successive scale heights in an arbitrary atmosphere, km: 

s h = S h 1 when P'co * 

Sh = Sh 2 when p M2 < P w < cT^ 

s h = s h 3 when p B > 


frequency factor for Arrhenius rate constant, K, for 
pyrolysis reaction (eq. (5)) 


t 


time, sec 



t G ff time, t, at which external loading is set to zero in computing 

programs; shut off time for wind-tunnel cases with calcu- 
lated loadings; for flight cases or arbitrarily programmed 
loadings, should be made larger than any time considered, 
sec 


TV 


ref 


A w 


wc 


wp 

T-i 


u 


u 


e 


v 


temperature, °K 

initial temperature deep in material, °K (eq. (122)) 

reference temperature, °K; may be used in material properties 
subprograms; generally zero 

wall (front face) temperature, °K 

corrected wall temperature., °K (see discussion following 
eq. (119)) 

predicted wall temperature, °K (eq. (119)) 

temperature array of preceding, time line based on the y_ \ 
array; used to calculate K _ \ array; °K 

horizontal component of trajectory velocity, km/sec 
(eqs. (88), (89)) 

tangential velocity at outer edge of boundary layer, cm/sec 
(eqs. (Bl) - (B3) ) 

vertical component of trajectory velocity (positive upward), 
km/sec (eqs. (88) -(90)); also normal velocity in boundary 
layer as used in appendix B, cm/sec 


v air 


v 


s 


v s comb> v sero’ 
v s sub 


V 

V 


velocity of external (ambient) gas normal to front surface 
and considered positive, cm/sec (eqs. (24), (25)) 

velocity of solid material in coordinate system (y direc- 
tion) , cm/sec; |v s | = rate of surface recession (eq. (2)) 

v s due respectively to combustion, erosion, sublimation; 
cm/sec (eq. (2)) 


volume occupied by pyrolysis gas per unit spatial volume, 
dimensionless (eqs. (15), (17)) 

enthalpy velocity, km/sec; defined as V 2 = 0.00836h st 
(eq. (56)) 

free-stream (trajectory) velocity, km/sec (eq. (88a)) 
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w rate of pyrolysis, from solid to gas per unit space volume, 

g/cm 3 sec (eqs. (3), (4)) 

x longitudinal coordinate along meridian (in boundary layer) , cm 

(eqs. (Bl) - (B3) ) 

y coordinate normal to front surface (inward) , cm 

y(0.10) depth of 10-percent char degradation (J = 0.90), cm 

y(0.90) depth of 90-percent char degradation (J = 0.10), cm 

y_i array of y values from the preceding time line, referring to 

the same points in the material as the current y array; 
deeper than the y array by the increase in surface reces- 
sion between the two time lines; y_i(M) = y(M) + AX, cm 
(see discussion following eq. (11)) 

z array defined by equation (9) 

Y_i array, same as z, but for the preceding time line and spaced 

according to the y_ 1 array (used in eq. (12)) 

Z x stability parameter for finite differencing, dimensionless; 

defined by equation (111); should be < 1/2. 


u 1 , CL 2 y 

U 3 y ^4 

y 

A 


A 


V 2 


e 


constants in equation (46) 


trajectory angle, deg; positive above horizontal (eq. (88b)) 

thermal thickness, defined by equation (96), cm; also increment, 
for example. Ah = enthalpy potential, cal/g 

increment of enthalpy potential due to homogeneous combustion 
(eq. (83)); in V 2 units; km 2 /sec 2 

surface emissivity, dimensionless (eq. (80)) 


6p activation temperature for the Arrhenius rate, K, for the pyrol- 

ysis reaction, °K (eq. (5)) 

y viscosity, poise (eqs. (B1)-(B7)) 

p density, g/cm 3 
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Pp,p c h> ) densities, g/cm 3 ; virgin plastic, pure char based on total 

> spatial volume (includes volume occupied by gas), pure char 

p cha ,p g >p s J based on volume occupied by solid particles of char only 

(eqs. (16), (17)), gas based on spatial volume (eq. (15)), 
solid (based on spatial volume) which can vary between p^ 
and p p (eqs. (7), (11)) 


P 


oo 


free-stream atmospheric density in Earth sea-level atmospheres, 
dimensionless , D/1226 


P 


CO 


2 


,p n 


3 


o 


values assigned to p^ at which changes of scale height in an 
arbitrary atmosphere occur; see Sft , s h 2 > s h 3 

Stefan constant, 1.369xl0“ 12 cal/cm 2 sec °K 4 (eq. (80)) 


4 > 


4 > 


approximation of the ratio of stored energy to the stored energy 
associated with an exponential temperature profile, 
dimensionless (eq. (97)) 

constant, dimensionless (defined by eq. (B7)) 


X 


surface recession, cm 


x l 


characteristic recession depth, cm, used in tumbling correction 
(eq. (61)) 


\jj convective heat blockage factor, dimensionless (eq. (78)) 

modified convective heat blockage factor, dimensionless 
(eq. (71)) 


Subscripts 

Subscripts listed below should be referred to when a subscripted symbol 
is not listed in the nomenclature above. Subscripts can be combined. 

air external gas or air 

BF back face 

c continuum; also corrected 

c' as defined by equation (113) 

calc calculated 

ch char 

comb combustion 
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d 


diffusion 


e outer edge of boundary layer; also quantity based on energy 

balance 

ero erosion (other than combustion or sublimation) 

ex excess 

FM free molecule 

g pyrolysis gas 

i initial 

L2,L3,L4 refer to the second, third, and fourth finite difference calcu- 
lated grid points; used in equations (117), (120); see 
appendix C, sketch (h) 

M,N index numbers 

max maximum 

MF in semi-infinite calculations, refers to greatest depth considered 

(same as BF) ; in finite depth calculations refers to a virtual 
point (whose depth may be greater than the BF point) ; see 
Numerical Analysis and Computation Procedures, Some Features of 
the Computer Program, Back Boundary Condition and appendix C, 
sketch (h) 

0 no blowing 

ox oxygen 

p virgin plastic; also predicted 

r reaction 

s refers to solid ablating material (which can vary from virgin 

plastic to pure char) 

st stagnation 

sub sublimation 

theo theoretical 

w 

1 
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wall (front face) 
dummy variable 


previous time line 
free stream 


Superscript 


time derivative 



APPENDIX B 


CALCULATION OF AIR FLOW FOR OXYGEN DIFFUSION 


The amount of oxygen that can diffuse to the solid surface to support 
heterogeneous combustion depends originally on the availability of oxygen 
which, in turn, depends on the rate of air flow into a control volume. We 
take the outer edge of the control volume to be at the assumed flame front for 
homogeneous combustion. The homogeneous combustion causes some depletion of 
oxygen, and it is necessary to calculate a corrected concentration of oxygen 
at this location which then (less the surface concentration) becomes the driv- 
ing potential for diffusion of oxygen to the surface. We use a Lewis analogy 
for the rate of oxygen diffusion, so we want the corrected concentration of 
oxygen to be evaluated where we evaluate the enthalpy for the enthalpy poten- 
tial, the "outer edge" of the boundary layer. 

From boundary- layer theory and an integration of the continuity equation, 
we obtain 


(Pv) e - 



(Bl) 


where v is considered positive in the direction of the surface, and F is 
the standard dimensionless stream function as defined in reference 12. Equa- 
tion (Bl) corresponds with equation (6.14) of reference 12 when there is no 
mass transfer at the surface. Equation (Bl) is rewritten 



(B2) 


A modified Newtonian approximation is used for the gradient of u e : 



(B3) 


where the constant C ^ allows 
ate p st with the perfect gas 
the viscosity as 


p + to be evaluated 

z 2 

law in terms of p t 


in atmospheres . We evalu- 
and T st and approximate 


„ 1/2 
Pst = c B T st 


(B4) 
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When we evaluate constants for air, we have: 


/ Cpp) w 

Cpv) e - Cpv) w - ( p a - V CBS) 

The second terms on both sides of equation (B5) should be equal. This is also 
indicated by equation (6.14) of reference 12, which is written for any posi- 
tion in the boundary layer. Accordingly, 



Equation (B6) can also be obtained directly from equation (6.14) of reference 
12. The value of F e is not overly sensitive to variations in F w . This is 
shown in an extended (unpublished) version of reference 13, and is also indi- 
cated in reference 14. Analysis of results for air presented in reference 12 

indicates that for a relatively cold wall, the product, y(pu) w / (pu) s -t F e> 


is rather insensitive to a variation in wall enthalpy, and for air it has a 
value of about 2.95. For convenience, we will define 


$ 



2.95 


(B7) 


Then equation (B6) becomes: 


m a ir 

calc 


( p air v air-U 


calc 



(24b) 


where we have changed the subscript on pv to e & ^ c in equation (24b) . In the 

absence of other_data, <j> is normally unity. For oxygen-containing gases 
other than air, <J> may be adjusted because the constant, 0.1566, in equations 
(B5) and (B6) applies to air. The evaluation of the external gas (air) flow 
obtained from equation (24b) is always compared with the evaluation from equa- 
tion (24a) , and the lesser value is used for the rate of air flow as given by 
equation (25) . 


The purpose of the calculation described is to obtain the oxygen concen- 
tration after some degree of depletion by homogeneous combustion. Clearly, 
there is some uncertainty in the calculation, because we are working with the 
"outer edge" of the boundary layer. In general, this uncertainty is not 
critical. We can rewrite equation (30) by combining into it equations (26), 
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(27) , and (28) and obtain 


C oxf " C ox 


m gw E 2 

™air 


(B8) 


For moderate rates of gas expulsion, m™, an error in m a ^ r can make a 
smaller relative error in C ox f. For larger rates of gas expulsion, the rela- 
tive error in C ox f need not be small, but C 0X £ will be small, and there 
will be little heterogeneous combustion anyway. As shown in equation (29), 
C 0X £ is assigned the value zero, when the apparent calculated value .from 
equation (B8) is negative. This occurs at the largest mg W values, and then 
an error in m a; : r is_ of no consequence for heterogeneous combustion. The 
value of unity for <j> in equation (24b) appears to be valid for ablation in 
air. It gives results that are consistent with data reported in figure 11 of 
reference 13. 
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APPENDIX C 


USE OF COMPUTING PROGRAMS 


The computing program is actually a group of programs for handling several 
types of problems. Two of the main programs accept arbitrary heating rates, 
one of these programs being for a semi-infinite depth material while the other 
is for a finite depth material (with a back boundary condition applied) . Two 
other main programs (finite and semi-infinite) use heating rates that are cal- 
culated in the programs themselves. A fifth main program is a continuation of 
the semi-infinite arbitrary heating-rate program which allows a calculation to 
be interrupted and resumed (with possibly a new grid spacing) . The continua- 
tion program need not be a continuation of any previous calculations. It 
allows one to read in arbitrary initial conditions including temperature and 
degradation profiles, Tj_(y) and J-[(y), and rate of gas generation, rhg± (y) . 

The main programs can be used for either wind-tunnel or flight calculations. 

Subprograms are used with each main program as indicated by the groupings 
below . 

Main Programs 

Main A Arbitrary Q, semi-infinite 

Main B Arbitrary Q, semi-infinite, continuation 

Main C Arbitrary Q, finite 

Main D Aerodynamic Q, semi-infinite 

Main E Aerodynamic Q, finite 

Starting Values Subprograms 

Subprogram SV1 Starting values, arbitrary Q 
Subprogram SV2 Starting values , aerodynamic Q 

Front Surface Subprograms 

Subprogram FSEB1 Front surface energy balance, arbitrary Q 
Subprogram FSEB2 Front surface energy balance, aerodynamic Q 

Back Surface Subprogram (Typical) 

Subprogram QBACK Back surface energy balance, finite heat sink 

Trajectory Subprogram 

Subprogram DER Trajectory equations, aerodynamic Q 

Printing Subprograms (Typical) 

Subprogram PRINTSI Print, semi-infinite 
Subprogram PRINTF Print, finite 
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Material Properties Subprograms (Typical) 

Subprogram ALPROPA Material properties, Apollo- type material, j = 0.5 

Subprogram ALPROPB Material properties, high-density phenolic nylon 

material, j = 0.5 

Subprogram ALPROPC Material properties, same as ALPROPB except 

j = read in (B13) 

The starting values and front surface subprograms are selected to go with 
a main program of either the arbitrary heating-rate type or the aerodynamic 
cally calculated heating-rate type. The back surface subprogram is used only 
with main programs for finite depth, C and E. The trajectory subprogram is 
used only for aerodynamic heating cases (main programs D and E) and is then 
only needed for flight cases. The printing subprograms go with semi-infinite 
calculations (main programs A, B, D) or finite calculations (main programs C 
and E) . The material properties subprogram is selected to correspond to the 
ablation material being used. For some new type of material, it may be neces- 
sary to write a new material properties subprogram, but this is easily done. 
Within the grouping of main and subprograms listed above, there are also 
library-type subroutines that perform such operations as differentiation, 
integration, and interpolation. 

In the program listing above, the back-surface subprogram, the printing 
subprograms, and the material properties subprograms should be considered as 
typical; an unlimited number of variations of these subprograms can be written. 
The back surface subprogram that is listed considers the backing material as a 
heat sink with an arbitrary finite heat capacity. The program solves 
equations (87b) and (87c) . 

Variations in the listed printing subprograms can be made to print out 
any new quantities desired. For example, by interpolation, one can print out 
temperatures at fixed positions in the material (say at thermocouple depths) . 
Also, the units of the printed quantities can be changed. These changes 
involve simple calculations to be performed in the printing subprogram, and 
they do not affect any of the calculations performed in the main and the other 
subprograms when new quantities (not in common) are used in the printing 
subprogram. 

In writing a new material properties subprogram, certain quantities must 
always be determined. The following material properties are generally func- 
tions of temperature, T, or they may be constant. These properties are listed 
as arrays, the index being the grid-point number, M, from 1 to MF + 1 . The 
temperature array of the previous time line T(M) is used to evaluate the 
temperature-dependent properties. These quantities are: Cp(CPP(M)), 

c p dTi (PINT (M) ) , c ch (CPCH(M)), f c ch dTi (CHINT(M) ) , c pg (CPG(M)), 

-T = dKr> = 

I C pg dT j (GINT (M) ) , Kp (CKDBP (M) ) , -Jf (DKDBP(M)), K ch (CKDBCH (M) ) , 
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dK ch 

j ■ (DKDBCH(M)), and h sr (HSCM)) (see appendix A). These quantities can be 

obtained by interpolation of tabular data either read in or put integrally 
into the beginning of the subprogram (and the tables need not be in common) . 

Or if analytical expressions are available for some of the quantities, they 
can be evaluated thusly by writing the appropriate equations into the subpro- 
gram. Two other arrays are determined in the subprogram, K(CK(M)) and 
TC x (CKM1 (M) ) . These can be obtained as functions of arrays of temperatures, 

T and T_i in any way desired, but they are generally calculated by an 
Arrhenius equation (5), with the constants being read in. 

A number of single-valued quantities are determined as functions of sur- 
face and external conditions. The quantities associated with heterogeneous 
combustion are * chcomb (FI) , equation (43) , m cbcomb (DIRML = VW) , equa- 

d max 

tion (44), m chcomb (RERML = P) , equation (45), and h scom t (HCS) , equa- 
re max 

tion (85). The quantity, h scomb , can be calculated as a function of surface 
temperature, T w , in this subprogram, but it is often assigned a constant 
read- in value (HCHREF) in the subprogram. The effective heat of combustion of 
pyrolysis gases, h ce (HCG) is evaluated in the subprogram according to equa- 
tion (81), and it thus depends on hg C . This latter quantity can, in princi- 
ple, be evaluated as a function of boundary- layer temperatures and the extent 
of oxidation, which will in turn depend on a number of variables. This is 
clearly not warranted in view of the simple flame front model assumed. The 
quantity, hg C , is most appropriately put into the subprogram as a read- in 
constant, HGREF. 

The rate of loss of surface material other than by combustion and sublima- 
tion (called erosion), m c h ero (F2 = TAUWP) , is evaluated in the subprogram. 

This may be in the form of an empirical expression such as equations (46) or 
some other evaluation. The energy associated with this process, h sero (CS) is 
usually zero (see eq. (85)). If h sero is a variable function of surface 
temperature, say, it should be evaluated in this subprogram. If h sero is a 
constant, it is read in as CS, and it need not be included in this subprogram. 

The loss rates of surface material due to sublimation are calculated in 
this subprogram. These quantities are ^ c h su b (FSUB = UDB(5)), equation (53), 
m c hsFM (FSUBRAT = UDB(4)), equation (48), and m chsd (FSUBDIF = UDB(6)), equa- 
tion (52). The energy of sublimation h sub (see eq. (86)) is usually approx- 
imated as a constant, in which case it need not be included in the subprogram, 
but it is read in as E20. If one wishes to make h sub a function of surface 
temperature, the evaluation should be performed in this subprogram, and for 
h sub , UDB(2) should be used (which is in common). 


Computing Program Options 

In addition to the major options determined by selection of a main pro- 
gram and accompanying subprograms, there are several options that can be 
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selected within the programs: 

1. Variation of density of the ablating material 

(a) The solid density, p s , at each depth, y, is not allowed to 
increase with time. This is only applicable to wind-tunnel con- 
ditions with constant heat loading. It is not necessary to use 
this option. KG = 1 . 

(b) The normal situation allows p s to vary freely. (This is valid 
for wind-tunnel cases also.) KG = 2. 

2. Calculation of internal degradation 

(a) The normal calculation follows the degradation of a point in the 
material . Each point of matter moves in the coordinate system 
(because of surface recession and because the origin is at the 
front surface). LL1 = 1. 

(b) An approximate calculation ignores the movement of solid material 
in the coordinate system. This speeds up computing time and is 
reasonable when the surface recession is negligibly small. 

LL1 = 0. 

3. Running conditions with aerodynamic heating (main program D or E) 

(a) Normal wind tunnel, KF = 1 . No rarefaction effects; uses 
equation (73) . 

(b) Rarefied wind tunnel, KF = 2. Takes account of rarefaction 
effects, equations (62)- (71) . (This includes all wind-tunnel 
cases, but the computing time is longer than with option (a).) 
Options (a) and (b) allow the wind tunnel to be shut off if 
desired; the calculations will continue while the model cools. 

(c) Entry flight with an initial temperature profile in the material 
that is computed as an exponential profile. With an assumed T w ^ 
(> T 0 ) , the initial values, and A^, are computed by the 
program. KF = 3 . 

(d) Flight with arbitrary initial values of atmospheric density, D-[ , 
and thermal thickness of an assumed exponential temperature pro- 
file in the body, . KF = 4. 

4. Planet and atmosphere for flight case with aerodynamic heating (main 

program D or E, with KF = 3 or 4) 

(a) Earth entry with ARDC atmosphere (approximated exponentially with 
three programmed values of scale height) . Earth radius pro- 
grammed at 6440 km, and g j = 980 cm/sec 2 . KC5 = 1. 

(b) Arbitrary planet (Rp^ and gp^) with exponential atmosphere having 
arbitrary scale height (initial, two intermediate, and final 
values). KC5 = 2. 


Damping Constants 

The damping calculation is given in equation (112) . When damping is not 
used, all damping constants should have the value unity. The read-in quanti- 
ties, KCT and NOP, are the time-line numbers (N) for the first and second 
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damping changes, respectively. The tabulation below lists the quantities 
damped and the first, second, and third damping coefficients. 


Quantity 

First damping 
coefficient 

Second 

Third 

T 

El 

B14 

CY 

ii»g(y, t) 

B6 

B7 

B6 

™chcomb 

B9 

BIO 

B9 

™chero 

Bll 

B12 

Bll 

*chsub 

Bll 

B12 

Bll 


PS 

CMO 

PS 


Nomenclature of Computing Program 

The nomenclature in the computing program is in symbolic FORTRAN 
language. Separate listings of input and output data are shown below. 


Input Data 

Input data are listed below in their order of card punching. Card for- 
mats are shown in sketch (g) . All input data are printed by the program in an 
initial readout (see Sample Case, below). A quantity listed as an option is 
defined in the section, Computing Program Options. 

Following the definition of a quantity, certain programs or options may 
be listed. For other program or option selections, the quantity is not needed 
or used, but the card formats must be maintained. 

The spacing sketch, time sketch, and printing sketch (sketches (h)-(j)) 
are referred to in the input listings below. The spacing sketch shows the 
y(M) spacing. The fine-spaced y increments are not used in the finite dif- 
ference solution of equation (1) . After the finite-differenced answers are 
obtained, the fine increments are used as interpolation points for a closer 
spacing near the front surface of the printed quantities (e.g., T(M)). See 
Output Data section. 
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I 



card 


4 

8 

12 

16 

20 

KF 

KG 

N 1 

N2 

NF 

KM3 

KOMI 

KCM2 

KCM3 

m 


32 

K4 

Ikon i 


KCMNF 


44| 
M 3l 
ILI 


48 1 
MF 
NOP 


56 


60 


J2 kC5 


J JARnAROn JJR2 
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68 


KCTKMIKM2 


JJCH 


72 


ALL NUMBERS IN 14 FORMAT (MUST BE RIGHT JUSTIFIED) 


A1 


C 5 
E5 
EI3 
CHI I 
DDT1 


DELTJ 


RHOP 


A5 


CMG 

CKDB 


COIN 

FI 


TE( l ) 


TE (97) 


Mil 


RJ (97)J 


CMRGVG( I ) 


CMRGVG(97) 


TTAR( I ) 


T TAR (297) 


HSTAR(I) 


10 

A2(0PEN) 

B6 

BI4 

C6 

E6 

£14 

VCINFI 

DDT2 

£19 

RHOCH 

A6 

PH I BAR] 
|CMCH(OPEN)| 
RIN 
F2 

TE(2 ) 

TE (98) 
RJ(2) 


RJ (98) 
CMRGVG( 2 ) 

CMRGVG(98) 

TTAR(2) 


HSTARC297) 


PT2AR( I ) 


PT2AR (2 97) PT2AR (298) PT2AR(299) 
OOCWARt I ) OOCWARl 2 ) QOCWAR( 3 ) 


| QOCWA RCT)| 
ORAR (!) 


QRAR (29 7) 
RAR(I) 


RAR (297) 


TTAR (298) 
HSTAR(2) 

|HSTAR(298)| 

PT2AR(2) 


19 

A3 

B7 

BI5 

C7 

E7 

EI5(0PEN)j 
GAMA I 
DDT3 
E20 

RHOCHA| 
CKI 
COW 
TOFF 
VCON 
FSUB 
TE(3 ) 

TE (99) 
RJ{ 3) 

RJ (99) 
CMRGVG(3)J 

CMRGVG(99) 

TTAR(3) 

TTAR (299) 
HSTAR (3) 

|HSTAR (299) 
PT2AR (3) 


ORAR (2) 

QRAR (298) 
RAR (2) 

RAR (298) I 


PT2AR(300) 
OOCWARt 4 ) 


bOCWAR(29G)OOCWAR( 299)1 


ORAR (3) 

ORAR (299) 
RAR (3) 

. 

RAR (299) 


28 

A4 

B8 

3I6(0PE: 

C8 

E8 

EI6 

TWI 

DYI 

E2I 

CS 

CK2 

CE 

E24 

QCS1N 

TT 

TE(4) 


RJ(4) 


CMRGVG( 4 ) 


TTAR(4 


HSTAR (4) 


PT2AR(4) 


ORAR (4) 

QRAR (300) 
RAR ( 4 ) 

RAR (300) I 


37 

46 

55 

64 72 

Bl 

B2 

B3 

B4 

B9 

BIO 

Bl 1 

BI2 

) Cl 

C2 

C3 

C4 

El 

E2 

E3 

E4 

E9 

EIO 

Ell 

EI2 

SIGMA 

RO 

OMCO 

OKBAR 

DEUW(OPEN 

ALLOW 

HST 

DCI 

RH02I 

OMO 

EI7 

EI8 

E22 

E23 

CJS 

CBAR 

ACT 

HSREF 

cx 

CY 

PS 

PT2 

CMBW 

CMO 

RG 

TREF 

HCHREF 

HGREF 

E25 

E26 

E27 

E28 

QCGIN 

QSUBIN 

STOR 

RESID 

PSI 




TE( 5 ) 

TE(6) 

TE(7 ) 

TE(8) 

RJ( 5) 

RJ(6) 

RJ(7 ) 

RJ(8) 

1 CMRGVGtS) 

CMRGVG(6) 

CMRGVGt 7 ) 

CKRGVG < 8 ) 

1 1 

TTAR(5) 

TTAR(6) 

TTAR(7) 

H 

> 

00 

HSTAR (5) 

HSTAR (6 ) j 

HSTAR (7) 

HSTAR (8) 

PT2AR (5) 

PT2AR (6) 

PT2AR ( 7 ) 

PT2AR (8) 

OOCWARt 5) 

00CWAR( 6 ) 

OOCWARt 7) 

QOCWAR ( 8 ) 

QRAR (5) 

QRAR (6) 

QRAR{7) 

QRAR (8) 

RAR (5) 

RAR (6) 

RAR (7) 

RAR (8) 






ALL NUMBERS IN E9 3 FORMAT (DECIMAL POINT NEEDED) 


2 

3 

4 

5 

6 

7 

8 

9 

10 
I I 
12 

13 

14 < 

15 
T I 

T 1 3 
Rl 

RI3 

Gl 

G 1 3 
A I 

A38 | 
Bl 

B38 

Cl 

C38 I 
Dl 

D38 | 
El 

E38 | 
FI 


F36 


Sketch (g).- 


Input card format. 


ALL 

PROGRAMS 


MAIN 

PROGRAM 

B 


MAIN 

PROGRAMS 
A, B, C 



Fine spaced 
e.g , LN= 4 


Ay, 


Grid point 
no. = M 1 


y L 4 


'iliUL 1 


(Ay 2 ) (Ay 3 ) 


e g. , K2 s 4 


'Ml M2 
(calculated) 


MF+ I 
(with 

finite depth 
calculation) 

, | 

I <Ay4 ’ i i I 

M3 MF-I MF 

(max = 99) 


eg , Ml s 17 


Ay 2 , Ay 3 , Ay 4 are increments used in finite difference equation 
A y[ = increment used for interpolation 


Sketch (h) • - Spacing. 


Nl, N2 are break points for changing At 


T ime line no = N 


| — At ( = DDT I 
I 

t = t j 


-At 2 = DDT 2-* Atj= DDT3— 


Nl 


N2 


NF 


At fJ At 2 , Af 3 are increments used in finite difference equation 


Sketch (i).~ Time. 


One time Ime 

KM I, KM2, KM3 are breok points for printing increments 
| (KCMI) (KCM2)j (KCM3)|( KCMF) | 

Grid point no I KMI KM2 KM3 MF(semi-co calculation) 

MF + I (finite calculation) 

KCMI, KCM2, KCM3, KCMF are grid-point printing increments 


Between time lines 

KNI, KN2 ore break points for printing increments 


(KCNI) 


(KCN2) I (KCNF) 


Time line no. I KNI KN2 NF 

KCNI, KCN2, KCNF are time line printing increments 


Sketch (j).- Printing. 


Card A 

(All numbers are integers in 14 FORMAT) 

KF Running condition option, main program D or E 

KG Variation of density of the ablating material option 

N1 Time-line number at which the finite time increment. At, changes from 

At} (DDT1) to At 2 (DDT2) (see sketch (i)) 

N2 Time-line number at which the finite time increment. At, changes from 
At 2 (DDT2) to At 3 (DDT3) (see sketch (i)) 

NF Final time-line number (see sketch (i)) 

K2 Defined by Ay 2 = (K2)Ayx. Can be > 1 (see sketch (h)); Ay 2 is the 

smallest Ay used in the finite difference equation; quantities in 
the Ay! spacing are obtained by interpolation from the finite- 
difference obtained values in the Ay 2 spacing. 

K3 Defined by Ay 3 = (K3)Ay 2 ; can be > 1 (see sketch (h)) 

K4 Defined by Ay 4 = (K 4 )Ay 3 ; can be > 1 (see sketch (h)) 

LN Increments of Ay 2 spacing over which Ay^ spacing exists; must be > 3 
(see sketch (h)) 

M2 Grid point at which space increment changes from Ay 2 to Ay 3 (see 
sketch (h)); (LN) (K2) + 1 = Ml < M2 < M3 

M3 Grid point at which space increment changes from Ay 3 to Ay 4 (see 
sketch (h)); M2 < M3 < MF 

MF Grid point at deepest point calculated. Maximum value = 99 (see 

sketch (h)). For semi-infinite calculations (main programs A, B, D) , 
this depth may be greater than actual material depth, and we require 
MF > M3 > M2 > Ml . For finite calculations (main programs C, E) , 
this is the initial depth of material; for these cases, we require 
initially that MF > M3 > M2 > Ml . Also, if K2 ^ 1, require that 
MF > Ml + 2 after shrinkage. 

J1 Order of interpolation for quantities, ~z _ 1 (M) and T_ \ (M) , for previous 
time line. Used in equation (12). In the computing programs, the 
value of z is advanced from one time line to the next. 

J2 Order of interpolation for calculated temperatures, T, at fine-spaced 
grid points (see sketch (h)) 

KC5 Planet and atmosphere option (main program D or E with KF = 3 or 4) 

KCT Time line, N, for first damping change 
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KM1 Grid point at which printing interval on one time line changes 

from KCM1 to KCM2 (see sketch (j)); 1 < KM1 < KM2 - KCM2 

KM2 Grid point at which printing interval on one time line changes 

from KCM2 to KCM3 (see sketch (j)); 

KM1 + KCM2 < KM2 < KM3 - KCM3 


KM3 


KCM1,KCM2 , / 
KCM3 , KCMF ( 

KN1 


KN2 

KCN1,KCN2, 

KCNF 

LL1 


Card B 

(All numbers are integers in 14 FORMAT) 

Grid point at which printing interval on one time line changes 
from KCM3 to KCMF (see sketch (j)); KM3 > KM2 + KCM3; also 
KM3 < MF - KCMF for main programs A, B, D; 

KM3 £ MF - KCMF - 1 for main programs C, E (using MF 
value after shrinkage) 

Printing intervals of grid points (see sketch (j)) (normally 
unity) 

Time-line number at which time-line printing interval changes 
from KCN1 to KCN2 (see sketch (j)) 

Time-line number at which time-line printing interval changes 
from KCN2 to KCNF (see sketch (j)) 

Time-line printing intervals (see sketch (j)) 


Internal degradation option; LL1 = 1, normally 


NOP Time line, N, for second damping change 

NAR Number in loading array; maximum = 300; main programs A, B, 

and C 


JJAR 

NAROP 

JJR2 


Initial order of interpolation in loading array; main pro- 
grams A, B, C 

Loading array printing option; NAROP = 1, prints array; = 0, 
does not print; main programs A, B, C 

Final order of interpolation in loading array; main programs 
A, B, C 


JJCH Number in loading array at which order of interpolation 

changes from initial (JJAR) to final (JJR2) ; main programs 
A, B, C 


All remaining cards have each number in E9.3 FORMAT, with eight (or 
fewer) numbers per card. 
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Card 1 


A1 Ai, equation (54); main programs A, B, C; also main programs D, E 
with KF = 3 or 4 

A2 (Open) 

A3 A3, equation (75) 

A4 A4, equation (62), main programs D, E; otherwise, open 

B1 iL if constant (and unevaluated in the material properties subprogram), 

otherwise open 

B2 if constant (and unevaluated in the material properties 

subprogram) , otherwise open 

B3 Sp, equation (5) 

B4 0p, equation (5) 

Card 2 

B5 Should be zero for all cases, except that B5 = t Q ££ can be used for 

constant loading cases (main programs A, B, C and D, E with KF = 1 
or 2) . B5 ^ 0 is not necessary for these cases, but when t < B5, 
it insures that K ^ K_i. 

B6 First and third damping coefficients, ihg(y, t) 

B7 Second damping coefficient, ihg(y, t) 

B8 T 0 

B9 First and third damping coefficients, n' c hcomb 

BIO Second damping coefficient, A cbcomb 

Bll First and third damping coefficients, ™ c hero anc * ™chsub 
B12 Second damping coefficient, m chero and m cbsub 

Card 3 

B13 j, equations (36)- (41), (45). ALPROPC; otherwise open 
B14 Second damping coefficient, T w 

B15 e 
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B16 (Open) 

Cl Ci, equation (60); main programs D, E 

C2 C 2 , equation (60); main programs D, E 

C3 C 3 , equations (92 ) 9 (95); main programs D, E with KF = 3 or 4 

C4 C 4 , equations (92), (95); main programs D, E with KF = 3 or 4 

Card 4 

C5 Sy^, equation (90); main programs D, E with KF = 3 or 4 and KC5 = 2 
C 6 Cg, equation (62); main programs D, E 

C7 L/Dr, equation (89); main programs D, E with KF = 3 or 4 

C 8 g i x 10" 5 , equation (89b); main programs D, E with KF = 3 or 4 and 

P KC5 = 2 

El First damping coefficient, T w 

E 2 E 2 , equation (27) 

E3 E 3 , equation (57b) 

E4 E 4 , equation (79); main programs D, E 

Card 5 

E5 Eg, equation (79); main programs D, E 

E 6 Eg, equation (79); main programs D, E 

E7 E 7 , equations (50), (51) 

E 8 Eg, equation (49) 

E9 Eg, equation (94); main programs D, E with KF = 3 or 4 

E10 E 1 0 , equation (57b) 

Ell M cv , equations (47), (48), (52) 

E12 E 1 2 , equation (91); main programs A, B, C; also main programs D, E 
with KF = 3 or 4 
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Card 6 


E13 

E14 

E15 

E16 

SIGMA 

RO 

OMCO 

OKBAR 

CHI 1 
VCINFI 

GAMA I 
TWI 

DELTW 

ALLOW 

HST 

DCI 

DDT1 

DDT2 
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Ex 3 , equations (93), (95); main programs D, E with KF = 3 or 4 
A C q, equation (64); main programs D, E with KF = 2, 3, or 4 
(Open) 

E x g , equation (61); main programs D, E with KF = 3 or 4 
o = 1 .369x10" 12 

R^, equation (60); main programs D, E 

(M/CdcA)., equation (95); main programs D, E with KF = 3 or 4 
K, equation (61); main programs D, E with KF = 3 or 4 


Card 7 


X l3 equation (61); main programs D, E with KF = 3 or 4 

V^.; main programs D, E; for flight cases (KF =3, 4); for wind- 
tunnel cases (KF = 1, 2), until t > t Q ££ 

y^; main programs D, E with KF = 3 or 4 

T w ^; all main programs, but with main program D or E and KF = 3, 
require T wi > T Q 

(Open) 

n, equation (4) 

h st ; main programs D, E. With KF = 1 or 2 (wind-tunnel case), h st 
is a constant until t > t Q ££, when h s t becomes zero; for the 
flight case with arbitrary initial conditions (KF = 4) , HST is 


D^; main programs D, E; with KF = 1 or 2 (wind-tunnel case), 

D = Di = constant until t > t Q ££ when D becomes zero; for the 
flight case with arbitrary initial conditions (KF = 4) , DCI is 
Di 


Card 8 


At x (see sketch (i)) 
At 2 (see sketch (i)) 



DDT3 

At 3 (see sketch (i)) 


DY1 

AYi (see sketch (h)) 


RH021 

Cpg if constant (and unevaluated in the material 
subprogram) 

properties 

OMO 

equation (74c) 


E17 

Epy, equation (60); main programs D, E 


E18 

Ei g, equations (92), (95); main programs D, E with KF = 3 or 4 


Card 9 


DELTJ 

R^l, equation (89); main programs D, E with KF = 
*KC5 = 2 

3 or 4 and 

E19 

Ei 9 , equation (49) 


E20 

h su bj equation ( 86 ) (positive quantity) 


E21 

A cv , equations (47), (48) 


E22 

c pwi’ see starting value equations (127), (134), 
programs D, E with KF = 1, 2, or 3 

(135) , (136) ; main 

E23 

Kpwi ; see starting value equations (123), (127), 
grams A, C and main programs D, E with KF = 1 

(136) ; main pro- 
, 2 , or 3 

CJS 

Cj, equation (14) 


CBAR 

c", equation (87c), main programs C, E 



Card 10 


RHOP 

Pp 


RHOCH 

Pch 


RHOCHA 

p cha 


CS 

h ser o» equation (85) 


ACT 

, equation (59); must be read in zero with main 

program E 

HSREF 

h sre f^ used to evaluate h sr in the material properties subprogram 
negative for an endothermic pyrolysis reaction at a reference 
temperature (usually 0 ° K) 
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CX c x , equations (43), (44), (45) 

CY Third damping coefficient, T w 


A5 

A 6 

CK1 

CK2 

PS 

PT2 

CMBW 

CMO 


Card 11 

A 5 , equation (77a) 

Converts pressure units from atmospheres to dynes/cm 2 = 1.013xl0 6 
Kj, equations (38), (39), (40) 

K 2 , equation (46) 

First and third damping coefficients, ip 


p t or p w , main programs D, E with KF = 1 or 2 (wind-tunnel cases) ; 


'Pt 2 or Pw 


is a constant for these cases 


Mp, equation (77a) 

Second damping coefficient, 


Card 12 


CMG Mg, equation (15) 

PHIBAR J, equations (24), (B7) 

COW C ox , equation (26) 

CE C E , equation (35) 

RG R , equation (15) (8.317xl0 7 ergs mole -1 °K -1 ) 

TREF T ref’ can usecl i n material properties subprogram (usually zero) 

HCHREF ^scomb if constant, equation (85); if h SCO mb is a function of 
temperature in the material properties subprogram, HCHREF is a 
reference value of hgcomb at some reference temperature 

HGREF hg C if constant, equation (81); if hg C is a function of tempera- 
ture in the material properties subprogram, HGREF is a reference 
value of hg C at some reference temperature 
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Card 13 


CKDB ^kdb 5 equation (43) 

CMCH (Open) 

TOFF t o£f’ w ^ en t - ^off* ^°° J h s t» and D become zero; shut-off time 

for wind tunnels (main programs A, B, C and D, E with KF = 1 or 
2); for flight cases (main programs A, B, C and D, E with KF = 3 
or 4) , t D ff must be > the flight time 

E24 "p^ ; main programs D, E with KF = 3 or 4 and KC5 = 2 

E25 "p ; main programs D, E with KF = 3 or 4 and KC5 = 2 

E26 , equation (90); main programs D, E with KF = 3 or 4 and 

2 KC5 = 2 

E27 , equation (90); main programs D, E with KF = 3 or 4 and 

3 KC5 =2 

E28 S^., equations (132), (133), and contained in equation (137); main 

1 programs D, E with KF = 3 and KC5 = 2; has arbitrary value but 
may equal (C5) 

Cards 14 and 15 and three sets of array cards listed below are used only 
with main program B (arbitrary Q, semi-infinite, continuation) . These addi- 
tional inputs are initial values that may be arbitrarily assigned or obtained 
from final values from main program A when a continuation calculation is being 
made. Cards 14 and 15 and the three sets of array cards (T^ (y) , J^y), and 
ihgi(y)) are all in the same format as cards 1 through 13. 

Card 14 


COIN 

Q w ^, equation (106a) 

RIN 

Qradi’ equation (106b) 

VCON 

Qvconi* equation (106g) 

QCSIN 

Q cs i> equation (106d) 

QCGIN 

Q C g^ , equation (106c) 

QSUBIN 

Qsubi’ equation (106e) 

STOR 

Qstori’ equation (106h) 

RES ID 

Q res i> equation (107) 
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Card 15 


FI 

™ch combi* equation 

(43) 

F2 

^cheroi* equation 

(46) 

FSUB 

™chsubi> equation 

(53) 

TT 

*i 


PSI 

equation (78) 



For a continuation calculation, it is also necessary to read in: 0M0 (Mp.), 

card 8, equation (74c) and ACT (Xj_) , card 10, equation (59). Following ckrd 
15, the initial temperature array, T^, is read in: TE (M) on cards Tl, T2, 

. . ., T13. The array of Ji is then read in: RJ (M) on cards Rl, R2, 

. . ., R13. Then the array of ihg is read in: CMRGVG (M) on cards Gl, G2, 

. . ., G13. The above three arrays have the dimension, MF, with a maximum 

value of 99 (up to 13 cards each) . The index M of each arrayed quantity 
corresponds to a depth y (array Y(M)). The spacing of the y array can be 

selected arbitrarily; it is not necessarily the same as that of prior calcu- 
lations (when a continuation calculation is being made) . 

Loading cards are required for main programs A, B, and C. For main pro- 
gram B, the loading cards follow the Gl, G2, . . ., G13 cards. For main pro- 
grams A and C, the loading cards follow card 13. The loading arrays have the 
dimension, NAR, with a maximum value of 300 (up to 38 cards each), and they 
are in the same format as cards 1 through 13. 

The loading arrays in order are: 

Time, t: TTAR (NN) on cards Al, A2, . . ., A38. 

Total enthalpy, h s -£ : HSTAR (NN) on cards Bl, B2, . . ., B38. 

Pressure, pp or p w : PT2AR (NN) on cards Cl, C2, . . ., C38. 

Cold-wall convective heating rate, q 0Cltf : QOCWAR (NN) on cards Dl, D2, 

. . ., D38. 

Radiative heating rate, q^: QRAR (NN) on cards El, E2, . . ., E38 . 

Effective nose radius, R: RAR (NN) on cards FI, F2, . . ., F38. 


Output Data 

The input data are printed in an initial output format (see Sample Case 
below). With main programs A, B, and C, the loading data arrays (t, h st , etc.) 
are printed or not when NAROP is assigned the value 1 or 0. The index num- 
ber of the loading arrays is also printed (labeled NN) . 

With main program B, the data read in on cards 14 and 15 are printed out 
with the first time-line output (always printed). (Some of these numbers may 
be slightly modified by being recalculated.) The quantities 0M0 (Mp^) and 

ACT (Xj[) are printed in the initial output, and they are also printed with the 
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first time-line outputs as TMLO and CHI, respectively. The arrays TE (M) , 

RJ (M) , and CMRGVG(M) are printed with the first time-line output. 

With all the main programs, the following additional quantities (some cal- 
culated) are printed in an initial output (see Sample Case below) . 

Q01I q w ^, equation (124) 

Q02I Tradi* ec L uat i° n (124) 

QOI q W ri > equation (124) 

TTI t^; equation (127) with main programs D or E and a wind-tunnel case 

(KF = 1 or 2) . Otherwise, t^ = 0 except with main program B (con- 
tinuation) when TT = t^ is read in. 

DELI Aj_; equation (123) except with main program D or E with KF = 4 

(when A^ is read in) . A^ is meaningful with an initial exponen- 
tial profile (eq. (122)), so it may not be meaningful with the 
continuation program, main program B. 

DCI Di; with main program D or E and KF = 1, 2, or 4, is read in; 

with KF = 3, is calculated from equation (137). With main 

program A, B, or C, is obtained from p t and V, 

equation (54) . 

RHO P ’00 = D/1226 (for this printing, D = D^) 

Ml Ml = 1 + (K2) (LN) ; see sketch (h) . 

The spacing of printing intervals is arbitrarily selected (see subsection 

Input Data and sketch (j)). For each time line printed out, the quantities 
listed below appear. (See subsection. Sample Case, below for format.) The 
first and last time lines (N = 1 and N = NF) are always printed out. 

N N, time-line number 

TT t 

R R 

SNGMA Sin y 

RHO F 

^ 00 

VC V 

RERML ™chcomb 

r max 

PSI ip 
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QOFM 

Hfm 

TMLO 

Mj 

QOC 

%c 

SQOO 


QO 

q 0 

PT2 

Pt 2 or Pw 

QOCW 

Q00 

^ocw 

Poo 

CHI+YDEG/CH 

X + y(0.90) 
original 

DIRML 

m chcomb 

YDEG/PL 

y(0.10) 


d max 



QR 

q r 

YDEG/CH 

y (0. 90) 

DTDYW 

0 T /9y)we 

PSIBAR 


OMC 

m/c d a 

C0IN 


TDEL 

A 

RIN 

Qrad 



B 

B 

CKTU 

K tu 

QSUB 

Psub 

STOR 

Qstor 

FSUB 

"'chsub 

VCON 

Qvcon 

FI 

m chcomb 

CHI+YDEG/PL 

X + y CO . 10) 
original 

F2 

m chero 

RESID 

Qres 

TWC 

T 

ERR 

e 

TWCAL 

T wp 

DTDYP 

(3T/ 3y) w 

XIP 

* = |v s i 


-l 

CHI 

X 

PSQO 

Qw 

DTT 

At 

FW 

^rad 

FSUBRAT 

A chsFM 

DTDYR 

OT/3y) wc 

FSUBDIF 

A chsd 

DSTOR 

Qstor 

PHI 


VARG 

Pvcon 

VC INF 


-VSSUB 

_v ssub 

ABROW 

lilrp 

DRES 

Pres 



-VS COMB 

_v scomb 


a depth from 
surface 


a depth from 
surface 
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-VSERO 


-v 


sero 


QCS 

q cs 


QCG 

q cg 


QCS IN 

Qcs 


QCGIN 

Qcg 


QSUBIN 

Qsub 


The four quantities below are printed out only by subprogram PRINTF (with 
main program C or E) as these are concerned with finite depth. 

QB 

^BF 


QBIN 

^BF 


DTDYBAK 

(9T/3y) B p 


TEMF 

T M p (a virtual point) 


Also for each time line, the 
in columns for the selected print 
Sample Case below for format) . 

following arrayed quantities are printed out 
spacing of grid points (see sketch (j) and 

M 

M (grid point number) 


Y 

y 


TE 

T 


RJ 

j 


RHOS 

Ps 


CMRGVG 

m g 


xz 

Z x , should be < 1/2 (for stability) 


For semi-infinite depth calculations (subprogram PRINTSI with main program A, 
B, or D) , the last arrayed quantities printed are for grid point number MF . 
For finite depth calculations (subprogram PRINTF with main program C or E), 
quantities for grid point number MF are not printed in the columns as MF is a 
virtual point. Quantities for grid point MF+1 are printed as this point 
represents the back surface of the material. 


79 


Sample Case 


The sample case used for illustration is a wind-tunnel case run under the 
conditions: total enthalpy, h st = 2700 cal/g; free-stream velocity, 

Voc = 3.0488 km/sec; pressure downstream of normal shock, p t ^ = 1.04 atm; for a 

duration of 2 seconds. The model has a hemisphere nose with a radius, R, of 
1.905 cm, and is Apollo type heat-shield material. This case could be handled 
by any of the five main programs; for example, with main program D or E, the 
loading would be calculated in the program. It was elected, for illustration, 
to use an arbitrary heating rate main program, so the loading was calculated 
independently and programmed in as an input. This loading is a cold-wall con- 
vective heating rate, q ocw = 226.2 cal cm” 2 sec" 1 and a radiative heating rate, 
qp 5 of zero. For further illustration, it was elected to use the continuation 
type main program B which uses a semi-infinite type calculation, and to start 
the calculation at a wind-tunnel time, t, of 1 second. The condition of the 
model at 1 second was previously obtained by use of main program A, and this 
condition (temperature and extent of degradation as functions of depth, rate 
of gas generation, extent of surface recession, etc.) furnished inputs to the 
continuation program. Calculations are for the stagnation point. The 
programs used for this example are: 

Main program B 

Starting values subprogram SV1 
Front surface subprogram FSEB1 
Printing subprogram PRINTSI 
Material properties subprogram ALPROPA 

The input data for the sample case with proper card formats are shown in 
figure 8. The print-out of the main part of the input data is given in 
figure 9. These are the data from cards A through 13 (common to all the main 
programs) . The loading data from cards A1 through FI are shown (common to 
main programs A, B, and C) , and the initial calculated values are also 
printed . 

The data read in from cards 14 and 15 and the arrays from cards T1-T5, 
R1-R5, and G1-G5 are peculiar to main program B. These data appear in the 
print-out of the first time line (N = 1), as shown in figure 10. All 
subsequent time-line print-outs are in the same format. 
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5 

10 15 

20 25 

30 35 

40 

45 50 

55 60 

65 70 

73 

CARD 

80 


2 4 1 5 1 

10 1 6 

1 1 4 

2 7 2 9 

3 9 1 

1 2 

1 2 

3 

A 


4 

1,1 1 

1 4 2 

5 8 10 10 

1 0 1 

5 1 4 

1 t 

1 1 0 


B 


1 .0 

0 . 0 

• 5 . 

o .o 

0 . 0 

0 . 0 

0 . 7 02 E 8 

1 10 60. 


1 


0 . 0 

1 .0 

! . 0 

3 11.0 

I . 0 

1 . 0 

1 .0 

1 . 0 


2 


0 . 0 

I .0 

. 7 5 

0 . 0 

0 . 0 

0 . 0 

0 . 0 

0 . 0 


3 


0 . 0 

0.0 

0 . 0 

0 .0 

1 . 0 

1.48 

.000084 

0 . 0 


4 


0 . 0 

0 . 0 

1 . 0 

2 2.2 9 3 

0 . 0 

0.1267 

3 3.0 

13.31 


5 


0 . 0 

0 . 0 

0 . 0 

0 . 0 

. 1 3 6 9 E - 1 

10.0 

0 . 0 

0 . 0 


6 


0 . 0 

0 . 0 

0 . 0 

3 11. 

0 . 0 

2 . 0 

0 . 0 

0 . 0 


7 


. 0 2 

. 0 2 

. 0 2 

.00 5 

0 . 0 

.0 1 8 8 0 2 

0 . 0 

0 . 0 


8 


0 . 0 

9 15 9 2. 

2 6 96. 

1 . 0 

0 . 0 

.0004547 

1 . 0 

0 . 0 


9 


.54 4 5 

.3 20 5 

2 . 25 

0 . 0 

.017014 

-83.3 

1 .3 3 3 3 

1 . 0 


1 0 


1.73 

1 .01 3 0 0 E 6 2 8 0000 . 

. 6 

1 . 0 


15.0 

1 . 0 


1 1 


1 5 . 

1 . 0 

. 2 3 

2 8 700 . 

8.3 1 700E70.0 

119 3. 

11 120. 


I 2 


2.05 

0 . 0 

2 000. 

0 . 0 

0 . 0 

0 . 0 

0 . 0 

0 . 0 


1 3 


14 6.8 

-19.273 

17.142 

6 .0 6 37 

2 7 .9 7 4 

- . 2 5 6 3 E - 

5 7 .6 8 8 3 

136.74 


1 4 


.010275 

' .000 5 5 06 5. 1 2 5 4 5 E - 

7 1.0 

.76012 





1 5 


2335.04 

1 2 10 7.11 

7 1879.194 

1651.270 

1423.346 

1195.422 

967.498 

90 1 . I II 


T 1 


834.723 

768 .3 36 

70 1 .9 4 9 

635.561 

56 9. 17 4 

545.318 

521 .461 

497.605 


T 2 


473.749 

449.893 

4 2 6 .03 6 

4 14 .588 

40 3.14 0 

391 .692 

300 . 244 

368 . 796 


T 3 


357.348 

327 . 357 

316.053 

312.375 

31 1 .332 

3 1 1 .072 

3 1 1 .014 

3 1 1 .002 


T 4 


3 11.000 

3 1 1 .000 

3 1 1.000 

3 1 1 .00 0 

3 1 1 .00 0 

3 1 1 .00 0 

311.000 



T 5 


. 27 0 5 6 E 

-4 . 5 3 9 9 4 E 

-4. I2582E- 

3.362I4E-3 

. ! 4 0 76 E - 

2 . 8 5 4 1 4 E - 

2.81 I82E-I0. 17320 


R 1 


O'. 3 5 9 3 9 

0.6 3 5 8 1 

0.8 67 3 4 

0 . 96824 

0 . 99375 

.99740 

0.99900 

0 . 99965 


R 2 ’ 


0.99989 

0 . 9 9 9 9 7 

0.99999 

1 . 0 

1 . 0 

i . 0 

i . 0 

i . 0 


R 3 


1 .0 

1 . 0 

1 . 0 

l'. 0 

1 . 0 

1 .0 

1 . 0 

1 . 0 


R 4 


1.0 

I . 0 

1 . 0 

1.0 

1 . 0 

1 . 0 

1 . 0 



R 5 * 


. 16 77 7 E 

- 1 . 1 6 80 8 E 

- i . 1 6 822 E - 

1 . 1 6 7 8 9 E - 1 

. 1 6 78 2 E - 

1 . 1 6 7 4 3 E - 

1 . 1 5 9 7 0 E - ( 

. 1 40 29 E - 

i 

G 1 


. I0 3 6 7 E 

- 1 . 5 5 3 5 5 E 

-2. I8775E- 

2’. 4 2 5 2 4 E - 3 

. 8 7 9 5 1 E - 

4 . 35 7 1 7 E - 

4 . 1 3 3 1 4 E -4 

. 4 62 2 7 E - 

5 

G 2 


. 1 4 3 3 8 E 

-50.0 

0 . 0 

0.0 

0 . 0 

0 . 0 

0 .0 

0 . 0 


G 3 


0 . 0 

0 . 0 

0 . 0 

0.0 

0 . 0 

0 . 0 

0 . 0 

0 . 0 


G 4 


0.0 

0 . 0 

0 . 0 

0 . 0 

0 . 0 

0 . 0 

0 . 0 



G 5 


o . o’ 

2 . 0 

2.001 

100.0 






A 1 


2 700 .0 

2 7 00 .0 

I 5 93.0 

15 9 3.0 






B 1 


1.04 

1.04 

0 . 0 

0 . 0 






C 1 


226.2 

2 26.2 

0 . 0 

0 . 0 






D 1 


0.0 

0 . 0 

0 . 0 

0 . 0 






El/ 


1.905 

1 .9 05 

1.905 

1.905 






F 1 

— 











— 


Figure 8.— Input data for sample case. 
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A 1 

A2 

A3 

A4 

B 1 


B2 

B3 


04 


B5 


Bt> 

07 

B6 

B9 

BIO 


011 

B12 


bl3 


B 14 


B2 5 

816 

Cl 

C2 

C3 


C4 

C5 


C6 


C7 


Ca 

El 

E2 

E3 

E4 


E5 

E6 


E7 


E 8 


E9 

E10 

Eli 

E 12 

El 3 


El 4 

E 1 5 


616 


SIGMA 


RO 

GMCO 

OK BAR 

CHI 1 

VC 1 NF I 


GAMA I 

T HI 


DELTW 


ALLOW 


HST 

DC I 

DOT 1 

DDT2 

DOT 3 


DY1 

RHC21 


OMO 


E 1 7 


El 8 

OELTJ 

E19 

£20 

621 


E22 

E23 


CJS 




RHGP 

RHGCH 

RHOCHA 

cs 

ACT 


HSREF 

CX 


CY 


A5 


A6 

CK1 

CK2 

PS 

PT2 


CMBM 

CMC 


CMG 


PHI0AR 


COW 

C£ 

RG 

TREF 

HCHREF 


HGREF 

CKDB 


CMCH 


TOFF 


0. lOOOGfc 

01 -o. 

0.53000E 03 

-0. 

-0. 

-0 


0.70200E 

08 

0.1 10606 

05 

-0. 


O.lOu^Ot 

..1 D.10000E 01 

U.311O0E 03 

O.IODDOE 01 

0.103006 n 

0 

10 JOOE 

31 0.100006 

01 

-0. 


0. 10000E 

01 

J.750u0t 

jO -j. 

'3 ■ 

0. 

3. 

0. 


-0. 


-0. 


0. 


-0. 

0.100Q0E 01 

O.140OOE 01 

0 . 84000E-04 

-0. 

-0 


-0. 


0.1OC00E 

01 

0.22293E 

02 

-u. 

0.12670E-00 

0.33003E 02 

0.1 331 OE 02 

-0. 

-0 


-0. 


-0. 


0.13690E- 

■11 

■J. 

-3. 

-0. 

-0. 

0. 

o. 


0.31 100 E 

03 

-0. 


0.200006 

01 

-0. 

-). 

3.23003E-01 

0.20000E-01 

0 .200006-01 

3< 

50000E 

-02 -0. 


0. 100O2E- 

-01 

0. 


0. 

-0. 

0.91592E 05 

0.2696CE 04 

0.100006 01 

-o 


0.4547CE- 

-03 

0.10000E 

01 



i). 54450E 

00 J.32050E-00 

0.22530E 01 

-0. 

0. 170146-01 

-0. 

83300E 

02 3.133336 

01 

0. luOOCE 

01 

3.17300E 

01 

0.101306 

<jl 3.28J0GE 06 

3 . 6 ,'OOCE 00 

O.IOOJOE 01 

-0. 

o. 

150006 

02 0. 13300 E 

01 

j. 15030E 

02 

O.IOOJOE 

01 

0. 230.0b- 

-33 3.287006 35 

0 . 33 173 E 00 

0. 

0.119306 04 

0, 

111206 

05 D.20500E 
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Figure 9.- Print-out of input data and initial calculated quantities for sample 

case. 



82 





N 

TT 

R 

SNGHA 

RHO 

VC 

RERHL 

PSI 


QOFM 

QOC 

GO 

QOCW 

GOG 

DIAHL 

OR 

DTDYW 

OHC 


TOEL 

CKTU 

QSUB 

FSUB 

FI 

F2 

TUC 

TWCAL 

XIP 


CHI 

DTT 

FSUBRAT 

FSUBDIF 

PHI 

VC INF 

ABROU 

TMLO 

SQOO 


PT 2 

CHl+YDfcG/CH 

YUEG/PL 

YDEG/CH 

PS IB AR 

COIN 

RIN 

B 

STOR 


VCON 

CHI +YDEG/PL 

RESID 

ERR 

OTDYP 

PSQO 

FW 

DTDYR 

DSTC3R 


VARG 

-VSSU6 

DRES 

-VSC0HB 

-VSERO 

QCS 

GCG 

QCS J N 

QCGIN 


QSUB I N 

1 

1.0000030 

0.1 9050 E 01 

0. 

0. 38014E-02 

0.47545E 01 

0. 16540E-OL 

0.7601ZE 

00 

0. 

0.16314b 03 

0.16314E 03 

0.22620E 03 

0. 

0 • 15375E-0 1 

0. 

-0 • 674 15E 05 

0. 


0.30024E-01 

0- 

-0.33821 E-04 

0. 12 545E-07 

0. 10275E-01 

0.5 5065E-03 

0.23350E 04 

0. 

0. 33777E- 

-01 

0. 170146-01 

0 . 200006-01 

0.47915E-O6 

D.12219E-07 

0.10030E 01 

0.30488E 01 

0.27603E-Q1 

0.18802E-01 

0. 


0.10400E 01 

0.40037E-J1 

0.5 16 18 E-D 1 

0.3 1023E-01 

0.76012E 00 

0 • 146806 03 

-0.19273E 02 

0. 34196E-00 

0.76883E 

01 

0.171426 02 

0.68632E-D1 

0. 13674E 03 

o. 

o. 

0.124006 03 

-0.30524E 02 

0. 

0. 


0.29745E 02 

J.39142E-U7 

0.55373E 02 

0. 320 59E-01 

0.17181E-02 

0.122586 02 

0.34752E 02 

0.60637E 01 

0.27974E 

02 

-0 .2 5630E-05 

M 

Y 

TE 

RJ 

RHOS 

CMRGVG 

KZ 




1 

0. 

2335.041 

0.27056E-04 

0.320516-00 

0. 167776-0 1 

0. 




2 

0.005 

2107.117 

0.539946-04 

0. 32051E-00 

0. 16808E-01 

0. 




3 

0.010 

1879.194 

0. 12 582E-03 

0. 320536-00 

0.1 682 2E-01 

o. 




4 

0.015 

1651.270 

0.36214E-03 

0.32058E-G0 

0.167 89E-01 

0. 




5 

0.020 

1423.346 

0.140766-02 

0.320826-00 

0. 16782E-01 

0. 




fa 

0.025 

1195.422 

0.8541 4E-02 

0 . 32241 E-00 

0.167436-01 

0. 




7 

0.030 

967.498 

0.811 82E-01 

C. 338686-00 

0. 15970E-01 

0 • 89637 E-Q 1 




8 

0.035 

901.111 

0. L732QE-00 

0.35930E-00 

3. 14029E-01 

0. 




9 

0.040 

834.723 

0.35939E-00 

0 .40 lOOE-OQ 

0.10367E-01 

0. 




10 

0.045 

768-336 

0 • 6358 IE 00 

0.46292E-00 

0.55355E-02 

0. 




11 

0.050 

701.949 

0.86734E 00 

0.51478E 00 

0. 18775E-02 

0. 




12 

0.055 

635. 56L 

0.96824E GO 

0.53739E 00 

0.42524E-03 

0. 




13 

0.060 

569.1 74 

0.99375E 00 

0 • 543 10E 00 

0.8795 1 E-04 

0.42557E-01 




14 

0.065 

545.3 18 

0 « 99740E 00 

0.54392E 00 

0.35717E-04 

0. 




15 

0.070 

521.461 

0.999D9E 00 

O.54420E 00 

0.13314E-04 

0 . 




16 

0.075 

497.605 

0.999656 00 

0.54442E 00 

0.46227E-05 

0 . 




17 

0.080 

473.749 

0.99989E OC 

0.544486 00 

0.143386-05 

0. 




18 

0.085 

449.893 

0.99997E 00 

0.54449E 00 

0 . 

Q. 




19 

0.090 

426.036 

0.99999E 00 

J.54450E 00 

J. 

0.51523E-01 




20 

0.095 

414.588 

O.IOOOOE 01 

0.54450E 00 

0 . 

0 . 




21 

0.100 

403.143 

0.10000E 01 

0.544506 00 

0 . 

0 . 




22 

0 • lu5 

391.692 

0.100306 01 

0. 544506 00 

0 . 

0 . 




23 

0.110 

380.244 

0.1003CE 01 

0.544506 00 

0 . 

0 . 




24 

0.115 

368.796 

o.ioooce oi 

0.54450E 00 

0 . 

0 . 




25 

0.120 

357.348 

O.IOOOCE 01 

0.5445CE CO 

0 . 

0 • 52402E-0 1 




26 

0.150 

327.357 

o.iooooe oi 

0.544506 00 

0 . 

0. 528046-01 




27 

0.180 

316.053 

o.iooooe oi 

0.544506 OC 

0 . 

0. 529576-01 




28 

0.210 

312.375 

o.iooooe oi 

0.544506 00 

0 . 

0.530076-01 




29 

0.240 

311.332 

o. icocoe oi 

0.5445CE 00 

0 . 

0.530216-01 




30 

0.270 

311 .072 

3.10000E 01 

0*544506 CO 

0 . 

0.5302IE-01 




31 

3.330 

311.014 

0. IQOJOE 01 

0.544506 00 

•j • 

0. 530216-01 




32 

0.330 

311.002 

o. looooe, oi 

0.544506 00 

3 . 

0.530216-01 




33 

0.360 

311.000 

o. looooe oi 

0.54450E 00 

0 . 

0 . 53C2 16-0 1 




34 

0.390 

311 .000 

0.100306 01 

0.544506 00 

0 . 

0.530216-01 




35 

0.420 

311 .000 

O.IOOJOE 01 

0.544506 00 

0 . 

0.530216-01 




36 

0.450 

311.003 

0.10CG9E 01 

J. 544506 00 

0. 

0.530216-01 




37 

0.480 

311.000 

o.iooooe oi 

0.544506 00 

c, . 

0.5302 1 E-0 1 




38 

0.510 

311 .000 

0.100006 01 

0.544506 00 

0 . 

0.530216-01 




39 

0.540 

311. OOj 

0.100036 01 

0.544506 30 

3 . 

u . 





Figure 10.— Print-out of first time line (typical) for sample case. 
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